Approved  for  public  release;  distribution  is  unlimited 


90  On 


Unclassified 


security  classification  of  this  page 


REPORT  DOCUMENTATION  PAGE 


1  a  Report  Security  Classification  Unclassified 

)b  Restrictive  Markings 

2a  Security  Classification  Authority 

3  Distribution  Availability  of  Report 

Approved  for  public  release;  distribution  is  unlimited. 

2b  Declassification  Downgrading  Schedule 

4  Performing  Organization  Report  Number(s) 

5  Monitoring  Organization  Report  Numberis) 

6a  Name  of  Performing  Organization 

Naval  Postgraduate  School 

6b  Office  Symbol 
(if  applicable)  33 

7a  Name  of  Monitoring  Organization 

Naval  Postgraduate  School 

6c  Address  (city,  state,  and  ZIP  code) 

Monterey,  CA  93943-5000 

7b  Address  (city,  state,  at.  i  ZIP  code) 

Monterey,  CA  93943-5000 

8a  Name  of  Funding  Sponsoring  Organization 

8b  Office  Symbol 
(if  applicable) 

9  Procurement  Instrument  Identification  Number 

8c  Address  (city,  state,  and  ZIP  code) 

10  Source  of  Funding  Numbers 

Program  Element  No  Project  No  Task  No  Work  Unit  Accession  No 

ll  Title  (include  security  classification'  MONTE  CARLO  GENERATION  OF  CERENKOV  RADIATION 


12  Persona!  Author(s)  Richard  J.  1 

Phillips  j 

13a  Type  of  Report 

13b  Time  Covered 

14  Date  of  Report  (year,  month,  day) 

1 5  Page  Count 

Master's  Thesis 

From  To  ■ 

December  1989 

122 

16  Supplementary  Notation  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official  policy  or  po- 
sition  of  the  Department  of  Defense  or  the  US.  Government. _ 


17  Cosan  Codes 

IS  Subject  Terms  (continue  on  reverse  if  necessary  and  Identify  by  block  number ) 

Field 

Subgroup 

Simulation  of  Gas  Cerenkov  Detectors 

19  Abstract  (continue  on  reverse  if  necessary  and  identify  by  blockjwmbw)- - —  — . .  * . ' 

the  Integrated  Tiger  Series  (ITS)  code  family, -produced  at  Sandia  National  Laboratories,  Albuquerque,  N.Mphnodel  the 
transport  of  electrons  photons.  These  codes  are  used  in  many  applications  such  as  radiation  shielding  and  radiation  dose 
prediction.  A  method  of  adding  Cerenkov  radiation  to  the  output  capabilities  of  the  code  was  devised  by  Joseph  Mack  of 
Los  .Alamos  National  Laboratory  and  Thomas  Jordan  of  Experimental  and  Mathematical  Consultants.  This  method  has 
been  extended  to  include  the  effects  of  wavelength  dependent  indicies  of  refraction.  The  capability  to  handle  air  and  carbon 
dioxide  gas  at  different  pressures  and  temperatures  has  been  added  to  the  program.  The  modifications  to  the  code  patch 
provide  the  user  with  wavelength  information  on  the  Cerenkov  spectrum  which  shows  44  percent  more  production  in  the 
1800  to  2000  angstrom  wavelength  bin  than  calculated  by  the  Mack-Jordan  patch.  These  additions  to  ITS  provide  a  poten¬ 
tially  valuable  tool  for  design  and  implementation  of  Cerenkov  detectors.  I'.  ,  i  '  r  l  .  //. , 

IT t  7  ■  i  o  ,  \  / 

,  /  ;  '  '  ■'  .  /  .  . 

/////>-  '  ' 

,  •  /  •  /  ■ 


20  Distribution  Availability  of  Abstract 

S  unclassified  unlimited  □  same  as  report  D  DT1C  users 

21  Abstract  Security  Classification 

Unclassified 

22a  Name  oi  Responsive  Individual 

X.K.  Maruyama 

22b  Telephone  t  include  Area  code  i 

( 408 1  646-2431 

22c  Office  Svmbol 

54  Ss 

DD  FORM  1473.84  MAR  83  APR  edition  may  be  used  until  exhausted  security  classification  of  this  page 

All  other  editions  are  obsolete 


Unclassified 


'  / 
/ 


9 


1 


Approved  for  public  release;  distribution  is  unlimited. 


Monte  Carlo  Generation  of  Cerenkov  Radiation 

by 

Richard  J.  Phillips 
Captain,  United  States  Army 
B.S.,  United  States  Military  Academy,  1979 

Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  PHYSICS 

from  the 

NAVAU  POSTGRADUATE  SCHOOL 
December  1989 


^  A 

^  Richard  J.  Phillips 


n 


ABSTRACT 


The  Integrated  Tiger  Series  (ITS)  code  family,  produced  at  Sandia  National  Labo¬ 
ratories,  Albuquerque,  NM,  model  the  transport  of  electrons/photons.  These  codes  are 
used  in  many  applications  such  as  radiation  shielding  and  radiation  dose  prediction.  A 
method  of  adding  Cerenkov  radiation  to  the  output  capabilities  of  the  code  was  devised 
by  Joseph  Mack  of  Los  Alamos  National  Laboratory  and  Thomas  Jordan  of  Exper¬ 
imental  and  Mathematical  Consultants.  This  method  has  been  extended  to  include  the 


effects  of  wavelength  dependent  indicies  of  refraction.  The  capability  to  handle  air  and 
carbon  dioxide  gas  at  different  pressures  and  temperatures  has  been  added  to  the  pro¬ 
gram.  The  modifications  to  the  code  patch  provide  the  user  with  wavelength  informa¬ 
tion  on  the  Cerenkov  spectrum  which  shows  44  percent  more  production  in  the  1800  to 
2000  angstrom  wavelength  bin  than  calculated  by  the  Mack-Jordan  patch.  These  addi¬ 
tions  to  ITS  provide  a  potentially  valuable  tool  for  design  and  implementation  of 


Cerenkov  detectors. 
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I.  INTRODUCTION 


Charged  particle  and  photon  transport  are  topics  of  much  interest  to  the  scientific 
community  and  have  been  for  much  of  this  century.  Since  Rutherford's  work  with 
scattering  from  thin  foils,  many  new  theories  have  been  developed  to  predict  the  radi¬ 
ation  and  charged  particle  deposition  caused  by  energetic  beams  of  charged  particles  or 
photons  diffusing  through  various  materials.  These  theories  describe  the  angular  de¬ 
flections  of  particles  and  the  energy  losses  associated  with  the  beam's  passage. 

Since  the  advent  of  the  high  speed  computer,  several  computer  codes  have  been  de¬ 
veloped  to  simulate  the  cascades  caused  by  energetic  electron  beams  impinging  on  dif¬ 
ferent  materials.  A  series  of  Monte  Carlo  codes  that  simulate  the  transport  have  been 
developed  by  Martin  Berger  and  Stephen  Seltzer  at  the  National  Bureau  of  Standards 
(NBS)  over  the  past  25  years.  These  codes,  such  as  ETRAN,  form  the  basis  for  other 
codes  written  at  laboratories  across  the  United  States.  The  ETRAN  algorithm  has  been 
used  by  Sandia  National  Laboratory  (SNLA),  to  develop  the  Integrated  Tiger  Series 
(ITS).  The  scheme  of  the  random  walk  used  in  the  ITS  calculations  has  its  roots  in 
ETRAN.  (Ref.  1] 

ITS  includes  S  different  variants  of  a  basic  code.  These  codes  reflect  various  com¬ 
binations  of  geometries  and  some  allow  the  superposition  of  electromagnetic  fields.  One 
of  these,  ACCEPT,  allows  the  user  to  specify  a  combinatorial,  fully  general,  three  di¬ 
mensional  geometry.  A  beam  of  electrons  or  photons  striking  and  diffusing  through  the 
target  is  simulated,  including  the  cascade  of  secondary  electrons  and  photons.  The 
program  tracks  primary  electrons  or  photons  and  all  possible  secondary  radiations,  in¬ 
cluding  knock-on  electrons,  bremsstrahlung.  Compton  electrons,  photoelectrons, 
electron-positron  pairs,  annihilation  radiation,  as  well  as  K-shell  characteristic  x-rays 
and  Auger  electrons.  Electrons  and  photons  with  energies  between  1  keV  and  1  GeV 
can  be  tracked.  [Ref.  1] 

Cerenkov  radiation  was  not  included  in  the  ITS  cascade  since  there  was  no  need  for 
its  inclusion  at  the  time  the  ITS  codes  were  developed.  However,  for  certain  applications 
such  as  the  design  of  Cerenkov  detectors,  its  inclusion  is  required.  This  problem  was 
addressed  by  Mack  and  Jordan.  They  produced  a  code  patch  to  allow  the  ACCEPT 
code  to  simulate  and  track  Cerenkov  production.  This  addition  to  the  code  samples  the 
Cerenkov  production  along  a  substep  and  then  ray  traces  (there  is  no  real  transport  for 
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these  low  energy  photons)  the  Cerenkov  light  through  the  geometry,  tallying  that  which 
arrives  at  the  detector,  and  the  Cerenkov  photons  produced  and  absorbed  in  each  zone. 
The  patch  due  to  Mack  and  Jordan  assumes  a  constant  index  of  refraction  for  each 
material  and  no  wavelength  dependence.  The  present  work  compares  the  results  of  this 
method  of  Cerenkov  sampling  with  the  results  produced  by  including  the  variable  index 
of  refraction  for  the  materials  in  which  the  Cerenkov  production  and  ray  trace  take 
place.  The  results  of  this  comparison  show  44  percent  more  production  in  the  1800  to 
2000  angstrom  wavelength  bin  than  calculated  with  the  constant  index  of  refraction  used 
for  runs  with  the  Mack-Jordan  patch. 
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II.  THEORY 


A.  HISTORY 

While  the  exact  date  of  the  discovery  of  Cerenkov  radiation  is  not  known,  we  do 
know  that  Mme.  Curie  was  among  those  to  notice  the  uncanny,  pale  blue  light  given 
off  from  bottles  of  concentrated  radium  solutions.  This  observation  was  documented  in 
1910.  In  1934,  Frank  and  Tamm  provided  the  correct  explanation  for  this  radiation. 
The  long  delay  between  the  time  when  the  radiation  was  first  noted  and  the  theoretical 
explanation  of  it  was  most  likely  due  to  the  abundance  of  new  discoveries  in  physics 
during  the  early  years  of  this  century,  especially  in  the  area  of  florescence  and 
phosphorescence.  Cerenkov  radiation  was  so  weak  as  to  cause  its  effects  to  be  masked 
by  the  presence  of  these  other  phenomena. 

Mallet  made  the  first  deliberate  attempt  to  study  the  phenomenon  in  the  years  1926 
through  1929.  However,  he  failed  to  pursue  the  subject  to  completion.  Then,  in  1934, 
Cerenkov  began  a  series  of  experiments  which  continued  into  1938.  His  experiments 
were  simple,  yet  they  showed  remarkable  agreement  with  the  theory  developed  by  Frank 
and  Tamm. 

In  1940,  Ginsburg  produced  a  quantum  theory  of  the  phenomenon  which  was 
thereafter  known  as  Cerenkov  radiation.  The  radiation  is  called  Cerenkov-Vavilov  in  the 
Soviet  literature.  In  1947,  Getting  proposed  using  a  photomultiplier  and  an  optical  sys¬ 
tem  for  detecting  single  particles  from  the  Cerenkov  radiation.  In  1951,  Jellev  was  the 
first  to  detect,  with  high  efficiency,  single,  fast  charged  particles.  His  work  was  followed 
closely  by  the  published  accounts  of  Mather  and  Marshall's  work  using  Cerenkov  de¬ 
tectors  for  the  measurement  of  velocities  of  particles  from  high  energy  accelerators. 

Since  that  time,  a  wide  variety  of  uses  have  been  found  for  Cerenkov  detectors,  in¬ 
cluding  the  discovery  of  the  anti-proton,  investigations  of  anisotropic  media  and 
ferromagnetics,  as  well  as  the  measurement  of  energies  released  in  nuclear  explosions. 

The  Cerenkov  effect  is  described  as  follows:  [Ref.  2:  pp.  5-6] 

1.  A  free  electron  injected  into  a  transparent  medium  produces  a  local  polarization 
of  the  medium. 

2.  At  small  velocities  this  polarization  is  roughly  spherically  symmetric.  This  local 
polarization  of  the  atoms  causes  them  to  behave  like  small  dipoles,  which  revert  to 
their  original  arrangement  after  the  passage  of  the  electron  (or  other  charged  par¬ 
ticle).  See  Figure  1  on  page  4.  However,  due  to  the  spherical  symmetry  of  the 
polarization  field,  there  is  no  resultant  field  at  large  distances. 
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Figure  1.  Local  Polarization  of  Transparent  Medium:  [Ref.  2:  pp.  5-6] 


3.  At  greater  velocities  of  the  electron,  the  spherical  symmetry  of  the  polarization  field 
is  destroyed  along  the  axis  and  radiation  is  emitted.  In  general,  this  radiation  in¬ 
terferes  with  itself  destructively  so  that  there  is  no  resultant  field  at  large  distances. 

4.  However,  if  the  velocity  of  the  electron  is  greater  than  the  phase  velocity  of  light 
within  the  medium,  then  there  is  an  angle  at  which  all  of  the  wavelets  from  all 
portions  of  the  electron  track  are  in  phase  with  one  another  so  that  there  is  a  re¬ 
sultant  field  at  a  distant  point.  See  Figure  2  on  page  5  and  Figure  3  on  page  6. 

The  angle  at  which  Cerenkov  radiation  is  produced  is  defined  by  the  Cerenkov  re¬ 
lation, 


(2.1) 


Cerenkov  radiation  is  only  produced  above  the  velocity  threshold  defined  by 


P 


min 


n  ’ 


(2.2) 


where  n  is  the  index  of  refraction  of  the  medium  and  /?  is  the  ratio  of  the  electron  velocity 
to  the  speed  of  light.  The  condition  for  Cerenkov  radiation  to  occur  is 


np>  1. 


(2-3) 


The  angle  of  emission  is  defined  bv 


6  =  cos'1  ( 

fin 


(2.4) 


Cerenko\  radiation  occurs  primarily  in  the  visible  and  near  ultraviolet  regions  of  the 
spectrum,  for  which  n  >  1  .  Emission  in  the  x-ray  region  of  the  spectrum  is  usually  not 
possible  since  the  index  of  refraction  is  then  less  than  unity  and  the  Cerenkov  condition 
cannot  be  satisfied.  (An  exception  is  at  absorption  edges,  where  it  is  possible  to  obtain 
n>  1  .) 

Cerenkov  radiation  may  be  thought  of  as  similar  to  the  V  shaped  bow  wave  from  a 
boat.  The  radiation  comes  off  at  a  particular  angle,  defined  by  the  index  of  refraction 
of  the  medium,  the  dominating  wavelengths  of  the  radiation,  and  the  velocity  of  the 
charged  particle.  See  Figure  3  on  page  6.  The  wavelength  of  the  radiation  plays  a  part 
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since  the  index  of  refraction  for  the  medium  is  dependant  upon  wavelength.  [Ref.  2:  pp. 
1-14] 

'  ’M  ■ 

B.  THE  CLASSICAL  THEORY  OF  FRANK  AND  TAMM 

In  1937,  Frank  and  Tamm  provided  the  first  theoretical  explanation  for  Cerenkov 
radiation.  Their  theory  made  the  following  assumptions:  [Ref.  2:  p.  15] 

1.  The  medium  is  considered  as  a  continuum,  so  that  the  microscopic  structure  is  ig¬ 
nored:  the  dielectric  constant  is  then  the  only  parameter  used  to  describe  the  be¬ 
havior  of  the  medium. 

2.  Dispersion  is  ignored,  at  least  in  the  first  approximation. 

3.  Radiation  reaction  is  neglected. 

4.  The  medium  is  assumed  to  be  a  perfect  isotropic  dielectric,  so  that  the  conductivity 
is  zero,  the  magnetic  permeability  equal  1,  and  there  is  no  absorption  of  radiation. 

5.  The  electron  is  assumed  to  move  at  constant  velocity;  ie.  the  slowing  down  due  to 
ionization,  and  the  multiple  Coulomb  scattering  are  ignored. 

6.  The  medium  is  unbounded  and  the  track  length  infinite. 
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Using  these  assumptions,  Frank  and  Tamm  derived  the  equation  for  the  output  of 
radiation  (cgs  units). 


1 

„2  2 
P  n 


(2.5) 


Without  a  cutoff  frequency  for  the  above  integral,  the  output  energy  is  infinite. 
However,  since  a  real  medium  is  always  dispersive,  the  radiation  is  restricted  to  bands 
such  that 


>  -J  ■ 


(2.6) 


Also,  the  absorption  bands  for  higher  frequencies  limit  the  radiation  to  the  near  ultra¬ 
violet  and  lower  frequencies.  The  refractive  index  is  less  than  unity  for  frequencies  in  the 
x-ray  region  making  Cerenkov  radiation  forbidden. 

For  constant  fin  we  can  express  the  radiation  intensity  in  terms  of  the  number  of 
photons  of  wavelengths  between  and  /2  which  are  emitted  by  an  electron. 


.V 


1  - 


(2.7) 


where  a  is  the  fine  structure  constant  and  /  is  the  path  length. 

There  are  a  variety  of  ways  to  express  the  spectral  distribution  of  Cerenkov  radi¬ 
ation.  The  intensity  of  light  of  frequency  v  may  be  written  as 

W=Xhv,  (2.8) 


where  A’  is  the  number  of  photons  of  that  frequency,  h  is  Plank's  constant  and  v  is  the 
frequency.  We  can  write  the  spectral  distribution  of  Cerenkov  radiation  in  the  following 
four  ways:  [Ref.  2:  pp.  15-22] 

•  Energy  per  unit  path  per  unit  frequency  interval. 


_£i r 

dldu) 
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•  Energy  per  unit  path  per  unit  wavelength  interval. 
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The  quantum  theory  of  Cerenkov  radiation  provides  more  detail,  but  in  accordance 
with  the  correspondence  principal,  yields  the  classical  results  for  macroscopic  phenom¬ 
ena.  Since  we  are  concerned  only  with  macroscopic  phenomena  in  this  work,  we  use  the 
classical  results.  The  constant  index  of  refraction  subroutine  uses  these  relationships  to 
calculate  the  amount  of  Cerenkov  radiation  produced.  Chapter  III  describes  this  ap¬ 
proach. 
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III.  THE  CONSTANT  INDEX  OF  REFRACTION  CERENKOV  PATCH 


Joseph  Mack  and  Thomas  Jordan  developed  a  set  of  sub-routines  for  the  ACCEPT 
[Ref.  3]  computer  code  which  give  the  code  the  ability  to  simulate  Cerenkov  production. 
Their  modification  calculates  the  number  of  Cerenkov  photons  A’  per  unit  path  length 
/  using  the  approximation, 


dS 

dl 


=  2  no. 


1  - 


nl  2 

P  n 


(>.2  >  >-i), 


(3.1) 


for 


/v 


dS 

dl 


2na 


h - 


1  - 


P2n\>.) 


(3-la) 


where  a  is  the  fine  structure  constant. 


Production  is  not  allowed  for  n  <  -~ 


A.  THE  Cl  PATCH 


The  ACCEPT  code  transports  electrons  using  multiple-interaction  pathlengths. 
This  means  that  the  code  simulates  the  electron  path  through  a  material  in  a  series  of 
substeps.  At  each  substep  the  electron  is  deflected  at  an  angle  which  is  sampled  ac¬ 
cording  to  the  multiple  scattering  theory  of  Gouldsmit  and  Saunderson.  (Ref.  4]  Energy 
losses  are  determined  from  the  Landau  distribution  [Ref.  5],  knock  on  electrons  are 
generated,  bremsstrahlung  photons  are  produced,  etc. 

The  constant  index  (Cl)  subroutine  links  Cerenkov  production  to  this  process.  It 
does  this  by  selecting  a  uniformly  random  point  along  the  pathlength  at  which  to  gen¬ 
erate  a  Cerenkov  photon  whose  weight  represents  the  expected  number  of  Cerenkov 
photons  produced  over  the  entire  pathlength.  ^See  Appendix  F  for  a  complete  dis¬ 
cussion  of  weight.)  The  emission  polar  angle  is  chosen  using  the  Cerenkov  relation, 
equation  (2-1).  The  azimuthal  angle  is  selected  using  a  uniform  random  distribution. 
The  mean  photon  energy  of  the  wavelength  band  is  calculated  and  assigned  to  the  sam¬ 
pled  photon. 
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Since  cross  sections  are  not  available  for  energies  as  low  as  those  associated  with  the 
Cerenkov  photons,  these  photons  are  not  transported  as  a  part  of  the  cascade.  Instead, 
each  Cerenkov  photon  is  ray  traced  through  the  geometry  using  geometrical  optics  and 
a  ray  tracing  algorithm  incorporated  into  the  subroutine.  These  ray  traces  depend  on 
the  geometry  and  on  a  variety  of  bulk  optical  properties  which  are  a  part  of  the  input 
file  supplied  by  the  user.  (Ref.  6]  These  inputs  are: 

1.  Index  of  refraction  for  each  zone.  This  number  is  the  wavelength  average  index 
of  refraction  for  each  material  in  the  problem. 

2.  The  absorption  coefficient  (cm-1)  is  used  to  calculate  absorption  losses  in  each  zone. 

3.  The  specular  reflection  fraction  is  the  fraction  of  light  that  undergoes  specular  re¬ 
flection  at  each  zone  boundary. 

4.  The  diffuse  reflection  fraction  is  the  fraction  of  light  that  undergoes  diffuse  re¬ 
flection  at  each  zone  boundary. 

5.  The  fractional  light  transmittance  is  the  fraction  of  light  transmitted  through  the 
zone. 

B.  CERENKOV  RAY  TRACE 

The  Cl  patch  works  as  follows.  The  user  defines  the  geometry  to  be  used  for  the 
problem  and  enters  this  as  a  list  of  zones  in  his  input  file.  (See  Appendix  E)  ACCEPT 
determines  the  length  of  a  substep  for  a  source  electron  as  it  enters  the  first  zone.  This 
is  a  pathlength  value  used  by  ACCEPT  to  calculate  interactions  with  the  material. 

If  the  Cerenkov  condition  is  satisfied,  equation  (3.1)  is  used  to  determine  the  total 
Cerenkov  weight  (Appendix  F)  produced  for  this  substep.  If  not,  the  routine  returns 
control  to  ACCEPT  for  the  next  electron  substep.  A  random  position  along  the  electron 
substep  is  sampled  and  used  as  the  starting  point  for  the  Cerenkov  ray.  The  Cerenkov 
angle  determined  from  equation  (2.1)  and  a  randomly  selected  azimuthal  angle  determine 
the  initial  path  of  the  ray. 

This  ray  is  traced  to  the  boundary  of  the  next  zone.  The  boundary  conditions  (the 
bulk  optical  properties  listed  above)  are  then  used  to  determine  whether  the  ray  is  diffuse 
or  specularly  reflected,  completely  or  partially  absorbed,  or  transmitted  into  the  next 
zone.  The  laws  of  geometrical  optics  and  Snell’s  law  are  used  to  determine  the  angles 
of  reflection  or  refraction  for  the  next  part  of  the  ray  trace. 

This  process  is  repeated  through  each  zone  until  a  termination  condition  is  reached. 
The  ray  can  only  terminate  in  three  ways: 

1.  It  gets  absorbed. 

2.  It  escapes  into  the  outside  universe  (the  last  zone). 
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3.  It  is  killed  because  it  gets  reflected  too  many  times.  This  parameter  (KALTOT)  is 
set  by  the  user  in  his  input  file.  For  the  problem  under  consideration,  KALTOT 
was  set  to  15.  This  kill  mechanism  is  included  to  eliminate  the  possibility  of  the  ray 
being  reflected  a  semi-infinite  number  of  times  between  surfaces  (geometric  insta¬ 
bility). 

Tallies  are  made  at  each  zone  to  record  the  Cerenkov  weight  produced,  the  energy,  the 
weight  absorbed,  the  weight  exiting,  and  the  weight  entering  the  zone.  The  code  also 
tallies  the  number  of  times  that  production  occurs.  These  tallies  are  saved  for  output 
and  of  course  one  of  the  zones  is  the  detector  zone  of  interest. 

This  model  is  quite  powerful  as  a  diagnostic  tool  for  analysis  of  Cerenkov  radiation. 
However,  it  assumes  that  the  indices  of  refraction  for  all  materials  are  constant.  In  the 
case  of  Cerenkov  production  in  a  gas  such  as  C02  the  index  of  refraction  is  very  nearly 
constant  over  the  entire  range  of  wavelengths.  For  example,  the  index  of  refraction  in 
C02  at  2  atmospheres  is  1.00110  at  1800  angstroms  [Ref.  7:  pp.432-434]  while  the  con¬ 
stant  index  used  in  the  Cl  routine  is  1.000894.  Despite  this  small  difference,  calculations 
using  a  table  of  index  of  refraction  by  wavelength  yield  44  percent  more  Cerenkov  pro¬ 
duction  in  the  1800  angstrom  bin  than  predicted  by  Cl.  The  patches  that  produce  these 
results  are  discussed  in  Chapter  IV.  A  detailed  comparison  of  the  results  of  these 
patches  and  the  reasons  for  the  differences  are  discussed  in  Chapter  V. 
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IV.  MODIFICATIONS  USING  VARIABLE  INDICES  OF  REFRACTION 


To  produce  a  code  patch  that  uses  the  variable  indices  of  refraction  for  the  material 
producing  the  Cerenkov  radiation  and  for  all  the  optics  through  which  these  photons 
pass,  two  different  approaches  are  taken.  They  are  the  multiple  trace  (MT)  method  and 
the  distribution  sampling  (DS)  method. 

A.  MULTIPLE  TRACE  METHOD 

The  MT  method  involves  splitting  the  spectrum  of  Cerenkov  photon  production 
into  small  wavelength  bins  across  which  the  index  of  refraction  (n)  is  very  nearly  con¬ 
stant.  The  wavelength  range  is  limited  to  1800-10000  angstroms.  This  range  is  subdi¬ 
vided  into  200-angstrom  bins  from  1800  to  5000  angstroms,  where  the  indices  of 
refraction  change  rapidly,  and  in  500-angstrom  bins  from  5000  to  10000  angstroms, 
where  the  indices  change  slowly.  A  plot  of  refractivitv  (n  —  1)  is  shown  in  Figure  4  on 
page  13.  This  shows  the  index  of  refraction  rising  rapidly  in  the  ultraviolet  region. 
However,  as  the  absorption  bands  arc  encountered,  the  index  of  refraction  falls  off 
sharply  until  the  Cerenkov  condition  is  no  longer  satisfied.  The  rapid  increase  in  values 
of  refractivitv  for  the  wavelengths  shorter  than  5000  angstroms  requires  the  use  of 
smaller  wavelength  bins  in  this  region.  Wider  bins  are  used  for  wavelengths  above  5000 
angstroms  to  optimize  code  efficiency. 

The  number  of  Cerenkov  photons  produced  in  each  wavelength  bin  per  pathlength 
dl  is  calculated  according  to  equation  (3.1).  Over  a  small  enough  range  of  wavelengths, 
the  index  of  refraction  is  approximately  constant  and  integration  of  (3.1)  over  wave¬ 
length  yields 


dS 

dl 


=  2rra(  1  — 


n2  2 
P  n 


L  A\ 


(/2 


(4.1) 


Equation  3.1  must  also  be  integrated  over  pathlength. 
homogeneous  material  this  becomes 


;Y  =  2jta/l 


r-ivAi 


'■2  J 


For  short  enough  substeps  in  a 

(>■2  (4-la) 


where  1  is  the  length  of  that  particular  substep.  (In  all  further  equations  where  ap¬ 
pears,  subsequent  integration  over  pathlength  occurs.) 
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Figure  4.  Refractivity  by  Wavelength:  This  figure  shows  the  index  of  refraction 
values  -  1  for  C02  in  the  MT  and  DS  methods. 
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The  number  of  Cerenkov  photons  per  wavelength  bin  is  converted  into  a  weight 
representative  of  all  photons  produced  along  the  entire  substep.  If  the  Cerenkov  con¬ 
dition  is  satisfied,  a  random  location  along  the  electron  substep  is  selected.  For  each 
wavelength  bin,  the  Cerenkov  angle  is  calculated  and  all  rays  of  appropriate  weight  are 
ray  traced  through  the  geometry  using  the  index  of  refraction  for  that  wavelength  and 
Snell's  law  in  each  of  the  optics.  Tallies  are  made  of  the  Cerenkov  weight  produced  in 
each  zone,  the  wavelength  of  the  light,  and  the  energy  produced.  The  Cerenkov  weight 
that  enters  and  exits  each  zone  is  also  recorded,  along  with  the  weight  that  is  absorbed. 
This  approach  is  illustrated  in  Figure  5  on  page  15.  This  shows  two  electron  substeps. 
Along  each  substep,  MT  samples  a  position  to  generate  all  of  the  Cerenkov  weight  for 
the  entire  substep.  It  then  calculates  the  weight  in  each  wavelength  bin  for  the  entire 
substep,  determines  the  Cerenkov  angle  for  each  wavelength  bin,  and  ray  traces  each 
weight  through  the  geometry.  (See  the  flow  diagram  for  the  MT  method  -  Appendix 
G.)  In  C02  the  range  of  Cerenkov  angles  when  is 

2.4°  <  0  <  2.7°.  (4.2) 

The  Cerenkov  weight  produced  from  a  substep  is  the  sum  of  the  weights  in  each 
wavelength  bin.  The  total  Cerenkov  weight  produced  is  then  the  sum  of  the  weights 
produced  integrating  over  each  substep.  If  the  total  of  the  weights  tallied  in  each 
wavelength  bin  were  different  than  the  sum  of  the  weights  produced  on  each  substep 
then  weight  balance  would  be  violated.  Clearly,  the  weight  balance  for  the  MT  method 
is  correct. 

The  output  from  the  code  patch  is  listed  by  wavelength  and  allows  the  user  to  an¬ 
alyze  the  spectrum  of  Cerenkov  radiation  reaching  each  zone.  The  update  file  which 
implements  this  method  is  included  as  Appendix  B. 

B.  DISTRIBUTION  SAMPLING  METHOD 

While  the  MT  method  provides  the  most  detailed  approach  to  examining  the 
Cerenkov  radiation  produced,  its  multiple  loops  cause  a  large  increase  in  computer  time 
over  the  Cl  method.  It  is  obviously  a  brute  force  method  that  forces  detailed  wavelength 
coverage.  With  this  in  mind,  a  method  which  samples  the  distribution  of  Cerenkov 
production  was  developed.  The  distribution  sampling  (DS)  method  is  intended  to  reduce 
the  amount  of  computer  time  required  to  obtain  a  similar  quality  result  to  that  provided 
by  the  MT  method. 


14 


1  0.000  angstroms 


28  Carankow  Rays 
from  samplad  position 
along  alactron  substap 
(Total  walght  producad  for 
this  substap  Is  dlstrlbutad 
among  tha  28  bins  in 
accordance  with 
aquation  4.1) 


Figure  5.  MT  Diagram:  This  figure  shows  Cerenkov  production  along  2  electron 
substeps  as  calculated  by  the  MT  method. 


Two  additional  subroutines  are  added  to  the  code  patch.  The  first,  CKVCDF.  cre¬ 
ates  a  discrete  cumulative  distribution  which  specifies  the  probability  for  Cerenkov  pro¬ 
duction  in  each  wavelength  bin.  This  routine  assumes  that  all  Cerenkov  production 
happens  in  the  interval  1800  to  10000  angstroms.  The  main  program  calls  CKVCDF 
once,  at  which  time  the  discrete  cumulative  distributions  for  all  the  materials  arc  calcu- 
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lated.  The  program  then  calculates  the  expected  number  of  Cerenkov  photons  produced 
over  the  entire  range  of  wavelengths  for  a  given  pathlength.  The  CDF  is  then  sampled 
through  a  call  to  the  routine  CKVWVL,  resulting  in  one  wavelength  bin  being  selected. 
The  total  expected  weight  is  assigned  to  that  sampled  wavelength.  This  Cerenkov  weight 
is  then  ray  traced  through  the  geometry  making  similar  tallies  as  in  the  MT  method. 

The  number  of  Cerenkov  photons  produced  in  each  wavelength  bin  is  unevenly 
distributed  as  can  be  seen  from  equation  (3.1).  This  is  shown  in  Figure  6  on  page  17. 
Each  step  on  the  CDF  plot  represents  the  cumulative  probability  that  the  sampled 
photon  is  in  one  of  the  bins  from  1800  angstroms  up  to  and  including  the  bin  represented 
by  the  step.  The  probability  that  the  photon  wavelength  is  sampled  in  any  particular 
bin  is  just  the  CDF  value  for  that  bin  less  the  CDF  value  for  the  previous  bin.  This  does 
not  mean  that  the  cutoff  for  Cerenkov  radiation  production  is  1800  angstroms.  This 
was  merely  the  cutoff  wavelength  selected  based  upon  available  index  of  refraction  data. 

More  photon  production  occurs  at  shorter  wavelengths.  This  calls  into  question  the 
validity  of  assigning  all  of  the  Cerenkov  weight  to  the  sampled  wavelength  bin.  How¬ 
ever,  this  disparity  in  production  is  accounted  for  in  the  discrete  cumulative  distribution. 
While  the  code  may  sample  a  wavelength  bin  in  the  longer  wavelength  region  and  ajsign 
too  much  weight  to  it,  the  probability  of  sampling  this  bin  is  smaller.  This  will  cause  a 
larger  sampling  frequency  at  shorter  wavelengths.  Providing  that  the  distribution  is 
sampled  enough  times  to  provide  valid  statistics,  this  approach  will  produce  results  sta¬ 
tistically  identical  to  those  of  the  multiple  trace  method,  while  utilizing  less  computer 
time. 

The  subroutine  CKVCDF  produces  the  cumulative  distribution  for  each  material 
as  follows.  The  total  number  of  Cerenkov  photons  produced  in  each  wavelength  interval 
(/  1  >  2 )  ? 


dl 


■  oczna 


(4.3) 


is  calculated  for  each  wavelength  bin.  These  individual  values  are  then  summed  to  yield 
a  normalizing  factor  (PTOT)  such  that 


PTOT  = 


(=i 


(4.4) 
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Figure  6.  Cerenkov  CDF:  This  figure  shows  the  discrete  Cerenkov  distribution 
sampled  by  the  DS  method. 

where  26  is  the  number  of  wavelength  bins  used. 
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The  discrete  probability  distribution  is  then  generated  by  letting 


)  _  '  / 
‘  PTOT 


(4.5) 


This  creates  a  normalized,  discrete  probability  distribution  (the  PDF)  with  F,  represent¬ 
ing  the  probability  of  generating  a  Cerenkov  photon  in  wavelength  bin  (i). 

To  form  the  discrete  cumulative  distribution  (CDF),  we  let 


26 

CDF{mtrlJ)  =  (4.6) 

/=  i 

where  mirl  is  the  material  number  and  /  is  the  index  for  the  wavelength  bin.  This  is  put 
into  a  loop  that  creates  the  CDF  for  each  material  and  stores  these  into  the  two  di¬ 
mensional  array  CDF(mtrl.i). 

This  CDF  array  is  randomly  sampled  (from  0  to  1)  to  determine  which  wavelength 
of  Cerenkov  light  will  represent  all  wavelengths  on  that  specific  electron  substep. 
Figure  7  on  page  19  shows  this  process  for  two  electron  substeps.  A  random  position 
is  chosen  along  the  substep,  all  of  the  Cerenkov  weight  generated  over  the  substep  length 
for  the  entire  range  of  wavelengths  is  assigned  to  a  photon,  and  the  wavelength  of  the 
photon  is  samp!  A'  from  the  discrete  cumulative  distribution  (CDF).  The  Cerenkov  angle 
is  calculated  for  this  wavelength  and  the  photon  is  ray  traced  through  the  geometry. 
The  next  electron  substep  is  taken  and  a  new  wavelength  is  sampled.  This  process  is 
repeated  for  each  substep.  (See  the  flow  diagram  for  the  DS  method  -  Appendix  H.) 

The  most  important  step  in  justifying  the  use  of  distribution  sampling  is  to  form  a 
weight  balance  equation.  We  must  show  that  the  sum  of  the  Cerenkov  weights  tallied 
in  each  wavelength  bin  equals  the  sum  of  the  Cerenkov  weights  produced  in  each  step. 
In  other  words  we  must  show  that  the  Cerenkov  weight  tallied  in  bin  (/)  is 

(4.7) 

where  w,  is  the  weight  produced  in  wavelength  bin  /  and  /  is  the  fraction  of  the  total 
weight  tallied  in  bin  (/).  IF  represents  all  the  Cerenkov  weight  produced  in  the  problem 
and  is  defined  by 
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Figure  7.  DS  Diagram:  This  figure  shows  Cerenkov  production  along  2  electron 
substeps  as  calculated  by  the  DS  method. 
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where  N  is  the  total  number  of  electrons  substeps  on  which  Cerenkov  production  occurs 
and  tpj  is  the  weight  produced  on  substep  j.  In  the  DS  method,  for  a  large  number  of 
substeps  f.Y)  ,  the  weight  produced  in  bin  i  is 


/V 

(4-9) 

7=i 

where  P,  is  the  discrete  probability  of  sampling  a  wavelength  in  bin  i  and  is  defined  by 
equation  (4.5 i.  As  the  number  of  substeps  on  which  sampling  occurs  (.V)  goes  to  infin¬ 
ity,  we  find 


and 


(4.101 


(4.11) 


Providing  that  S  is  sufficiently  large  to  provide  the  required  degree  of  confidence,  we  are 
justified  in  using  the  approximation  represented  by  equation  (4.11).  Substituting 
equation  (4.S)  into  (4.11)  yields  equation  (4.7).  This  shows  that  sampling  the  distrib¬ 
ution  provides  the  required  weight  balance  in  each  wavelength  bin  provided  that  the 
number  of  samplings.  N,  is  large  enough.! 

The  update  file  required  to  implement  this  routine  is  included  as  Appendix  C. 


!  An  excellent  discussion  of  the  estimation  of  Monte  Carlo  accuracy  is  given  in  the  MCNP 
manual.  | Ref.  8:  pp.  102-1 1 4 j 
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V.  RESULTS 


Runs  were  made  using  the  Cl,  MT,  and  DS  methods  to  model  a  gas  Cerenkov  de¬ 
tector.  The  geometry  for  this  detector  is  shown  in  Figure  8. 
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Figure  8.  Cerenkov  Detector  Geometry.:  This  figure  shows  the  geometry  for  a 
gas  Cerenkov  detector  that  was  used  to  compare  the  Cl.  N1T.  and  DS 
methods. 
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All  of  the  code  patches  calculate  the  total  Cerenkov  radiation  produced  in  each 
zone,  the  amount  entering  or  leaving  the  zone,  and  the  amount  absorbed  by  the  zone. 
However,  the  two  modified  patches  also  yield  data  on  the  wavelengths  of  radiation  in 
each  zone.  For  a  discussion  of  zones  and  geometry,  see  Appendix  E. 

A.  VALIDATION 

The  first  step  in  comparing  the  methods  was  to  insure  that  the  code  with  the  new 
patches  produced  similar  results  using  the  same  indicies  of  refraction.  To  check  this,  the 
variable  indicies  of  refraction  in  the  MT  method  were  replaced  with  the  Cl  constant 
indicies.  The  Cerenkov  spectrum,  as  calculated  by  the  Cl  and  the  MT  methods,  is 
shown  in  Figure  9  on  page  23.  The  output  of  the  Cl  and  MT  methods  are  identical. 
These  results  verify  that  the  Cerenkov  production  using  the  MT  method  reverts  to  the 
constant  index  (Cl)  method  when  the  same  indicies  of  refraction  are  used. 


B.  COMPARISON  OF  CODE  PATCHES 

The  first  area  of  interest  in  comparing  the  different  methods  is  the  amount  of 
Cerenkov  radiation  produced  in  zone  5,  the  C02  gas.  (See  Appendix  E  for  zone  de¬ 
scriptions).  Figure  10  on  page  24  shows  a  comparison  of  the  production  calculated  by 
the  constant  index  and  multiple  trace  routines  when  the  index  of  refraction  is  allowed 
to  vary  for  the  MT  method.  The  MT  method  produces  44  percent  more  Cerenkov 
weight  in  the  1800-angstrom  bin  than  the  Cl  method.  (For  this  problem  the  index  of 
refraction  used  for  Cl  in  each  of  the  materials  was  the  index  for  light  of  5000  angstroms.) 
The  MT  method  also  shows  about  5  percent  less  production  for  the  long  wavelengths 
than  calculated  by  Cl.  Figure  11  on  page  25  shows  the  same  comparison  as  a  ratio  of 
the  amounts  calculated  by  each  method.  The  uncertainties  in  the  the  ratio  plot  are  cal¬ 
culated  using  [Ref.  9:  pp.  56-89] 


2 


r 


2 


(5.0) 


where  r  =  and  a 2  is  the  variance.  The  ratio  plot  emphasizes  the  differences  in  the 
behaviors  oi  the  two  models.  Where  the  production  by  the  two  methods  is  the  same  on 
the  graph  is  where  the  MT  variable  index  for  C02  equals  the  constant  index  used  in  Cl. 
This  figure  illustrates  the  fact  that  when  the  Cerenkov  detector  photo-multiplier  tube 
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Figure  9.  MT  vs.  Cl  Cerenkov  Weights:  This  figure  shows  the  Cerenkov  weight 
produced  in  the  C02  as  calculated  by  the  MT  and  Cl  methods  when  the 
same  indicies  of  refraction  are  used. 
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Figure  1 1.  \1T  vs.  Cl  Production  Ratios:  This  figure  shows  the  ratio  of  Cerenkov 

weights  produced  in  C02  as  calculated  by  MT  to  that  of  Cl. 


has  a  non-constant  wavelength  sensitivity,  the  inclusion  of  the  proper  wavelength  de¬ 
pendent  index  of  refraction  is  important. 

Figure  12  on  page  27  shows  the  Cerenkov  production  calculated  by  the  DS  and  Cl 
methods.  Figure  13  on  page  28  shows  the  ratio  comparison  of  these  two  methods.  The 
same  qualitative  result  is  obtained  for  the  DS  and  MT  methods.  This  shows  the  larger 
production  at  the  short  wavelengths  and  smaller  production  at  the  long  wavelengths, 
although  there  are  large  statistical  variations  in  the  sampling  at  the  wavelengths  longer 
than  5000  angstroms.  The  degree  of  the  variation  depends  on  the  number  of  batches 
used  in  the  run  or  on  any  biasing  techniques  applied.  Biasing  has  not  been  used  to  im¬ 
prove  the  statistical  accuracy  for  this  problem. 

A  comparison  of  the  ratio  of  the  distribution  sampling  (DS)  to  the  multiple  trace 
(MT)  calculation.  Figure  14  on  page  29,  shows  further  evidence  of  this  statistical  fluc¬ 
tuation.  Figure  14  on  page  29  also  shows  that  the  results  of  DS  and  MT  are  nearly 
identical.  This  is  not  surprising  as  the  DS  method  is  designed  using  the  same  physics, 
but  incorporating  a  Monte  Carlo  sampling  technique  in  order  to  trace  fewer  Cerenkov 
rays,  thereby  improving  efficiency. 

The  next  important  zone  in  which  to  compare  the  methods  is  zone  18.  This  zone 
represents  the  photo-multiplier  or  detector.  An  analysis  of  the  Cerenkov  spectrum  here 
allows  the  user  to  see  what  Cerenkov  radiation  will  arrive  at  the  detector.  Figure  15  on 
page  30  shows  a  plot  of  the  Cerenkov  weights  reaching  the  detector,  as  calculated  by 
each  of  the  methods.  Note  that  the  Cl  method  only  calculates  a  total  weight.  The 
fraction.  of  the  total  weight  in  each  wavelength  bin  is  then  calculated  using 


(±jt)  l 

—  - 

1800  10000  ) 


(5.1) 


Figure  16  on  page  31  plots  the  ratio  of  the  Cerenkov  distribution  as  calculated  by  the 
MT  method  to  that  calculated  by  Cl.  A  plot  of  the  ratio  of  DS  to  Cl  weights  is  shown 
in  Figure  17  on  page  32.  Again  the  statistical  errors  in  the  DS  results  are  larger  than 
10  percent  for  the  wavelengths  longer  than  5000  angstroms,  although  DS  closely 
matches  the  MT  results. 

In  order  to  determine  whether  the  spectrum  arriving  at  the  detector  is  different  than 
that  which  is  produced,  it  is  necessary  to  examine  the  ratio  of  weight  that  arrives  at  the 
detector  to  weight  produced  in  the  C03.  Figure  18  on  page  33  shows  this  ratio  for  the 
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Figure  12.  DS  vs  Cl  Cerenkov  Production:  This  figure  shows  the  Cerenkov 
weight  produced  in  C02  as  calculated  by  DS  and  Cl. 
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Figure  13.  Cerenkov  Production  Ratio,  DS  to  Cl:  This  figure  shows  the  ratio  of 
the  Cerenkov  weight  produced  as  calculated  by  the  DS  method  to  that 
of  Cl. 
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Figure  14.  Cerenkov  Prodcution  Ratio,  DS  vs.  MT:  This  figure  shows  the  Ratio 
of  the  Cerenkov  weight  produced  in  the  C02  as  calculated  by  DS  to  that 
calculate  by  MT. 
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Figure  15.  Cerenkov  Weights  at  Detector:  This  Figure  shows  the  Cerenkov 
weight  reaching  the  detector  (zone  18)  in  each  wavelength  bin  as  cal¬ 


culated  by  the  MT.  DS.  and  Cl  methods. 
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Figure  16.  Ratio  at  Detector,  MT  to  Cf :  This  figure  shows  a  comparison  of  the 
Cerenkov  spectrum  reaching  the  detector  (zone  18)  as  calculated  by 
\1T  and  Cl. 
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Figure  18.  Ratio  Reaching  Detector,  MT  Method:  This  figure  shows  the  ratio  of 
the  Cerenkov  Weight  arriving  at  the  detector  to  that  produced  in  the 
CO  as  calculated  by  the  MT  method. 
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MT  method.  This  plot  shows  that  for  this  pan'eular  Cerenkov  detector  geometry,  the 
Cerenkov  wavelength  distribution  at  the  detector  is  roughly  the  same  as  that  which  is 
produced,  although  only  about  37  percent  of  the  photons  produced  make  it  to  the  de¬ 
tector.  This  is  because  some  of  the  radiation  is  absorbed  in  the  walls,  some  is  killed  from 
too  many  reflections,  and  some  escapes. 

The  ratio  of  the  amount  produced  to  the  amount  that  reaches  the  detector  for  the 
DS  method  is  shown  in  Figure  19  on  page  35.  Again  only  37  percent  of  the  photons 
produced  arrive  at  the  detector. 

The  Cl  method  assumes  that  the  distribution  arriving  at  the  detector  is  the  same  as 
the  original  distribution  and  it  shows  statistically  the  same  percentage  of  the  original 
radiation  reaching  the  detector,  35  percent. 

Finally,  Figure  20  on  page  36  provides  a  zone  by  zone  comparison  of  the  total 
Cerenkov  radiation  entering  the  zone,  as  calculated  by  each  of  the  methods.  While  the 
geometry  of  the  detector  is  a  complicated  set  of  19  zones  (described  in  Appendix  E),  the 
most  important  information  in  this  figure  is  the  Cerenkov  weight  arriving  at  the  detector, 
Zone  18.  Some  of  the  zones  show  more  Cerenkov  radiation  entering  than  was  produced. 
This  is  because  reflections  may  cause  Cerenkov  rays  to  enter  a  zone  more  than  once. 
The  total  Cerenkov  weight  is  conserved  if  the  Cerenkov  weight  produced  equals  the  sum 
of  the  weight  absorbed,  the  weight  killed  due  to  too  many  reflections  (more  than  the  user 
specified  number  of  allowed  reflections),  and  the  weight  that  escapes  into  the  outside 
universe  (the  last  zone.  19).  The  DS  and  MT  patches  both  include  routines  to  verify  that 
the  weight  conservation  fraction  is  unity.  Figure  20  on  page  36  provides  another  illus¬ 
tration  that  the  physics  used  in  each  of  the  codes  is  the  same,  although  the  methods  of 
calculation  are  different. 

C.  DEPENDENCE  ON  INDEX  OF  REFRACTION 

The  fact  that  there  is  44  percent  more  Cerenkov  radiation  produced  in  the  carbon 
dioxide  gas  as  calculated  by  MT  than  that  calculated  by  Cl  is  somewhat  surprising.  The 
index  of  refraction  for  carbon  dioxide  gas  at  2  atmospheres  is  1.00110  for  the  1800 
angstrom  bin.  The  constant  index  of  refraction  used  by  Cl  is  1.000894.  The  maximum 
difference  between  indicies  of  refraction  used  by  the  Cl  routine  and  the  variable  index 
routines  amounts  to  only  0.000206,  or  2  parts  in  10000. 

To  understand  how  such  a  small  change  in  the  index  of  refraction  can  produce  more 
than  a  forty  percent  increase  in  the  amount  of  Cerenkov  radiation  calculated  for  the 
shortest  wavelengths,  it  is  necessary  to  analyze  the  Cerenkov  production  formula,  l  or 
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Figure  19.  Ratio  Reaching  Detector,  DS  Method:  This  figure  shows  the  ratio  of 
the  weight  arriving  at  the  detector  to  that  produced  as  calculated  by  the 
DS  method. 
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Figure  20.  Total  Cerenkov  Weights  Entering  Each  Zone:  This  figure  shows  the 
total  weights  calculate  by  each  method. 
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a  small  enough  wavelength  bin  and  constant  n,  the  number  of  photons  produced  per 
path  length  is 


d.X 

dl 


=  2nca 


(5.2) 


For  the  extremely  relativistic  case  ,  /tel,  the  number  of  photons  per  path  length 
produced  in  this  bin  is 


dX 

dl 


(5.3) 


If  we  define  n  =  l  +  y,  then  this  becomes 


dX 

dl 


oc 


1  -• 
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(l  +  xY 


(5.4) 


In  the  case  of  carbon  dioxide  gas,  *<?1.  Therefore,  taking  the  first  order  Taylor  expan¬ 
sion  is  justified  and  the  above  proportionality  becomes 

~~  oc[i  —  (I  —  2x)].  (5.5) 


or 


dX 

dl 


zz2nx(2y) 


(5.6) 


Thus  a  change  from  y  —  0.000894  to  y  =  0.001 10  amounts  to  a  23  percent  increase 
and  it  produces  a  corresponding  23  percent  increase  in  the  amount  of  Cerenkov  radi¬ 
ation  calculated  for  the  1800-angstrom  bin.  Both  the  MT  and  DS  methods  calculated 
a  44  percent  increase  in  the  amount  of  Cerenkov  radiation  for  this  bin.  This  difference 
in  production  is  accounted  for  by  the  fact  that  the  constant  index  of  refraction  used  in 
the  Cl  routine  for  this  problem  is  below  the  average  index  of  refraction.  The  effect  of 
the  constant  index  being  below  the  average  index  can  be  removed  by  normalizing  the 
weights  produced  in  both  methods.  This  is  done  by  dividing  the  weight  in  each  bin  by 
the  total  weight  calculated  for  that  zone.  Figure  21  on  page  38  shows  that  once  this 
normalization  is  taken  into  account,  the  ratio  of  the  Cerenkov  weights  produced  is  about 
1.23  at  1800  angstroms. 
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Figure  21.  Normalized  Cerenkov  Production:  This  figure  shows  the  ratio  of  the 
normalized  Cerenkov  weights  produced  by  MT  to  that  of  Cl. 
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Another  important  characteristic  of  the  parameter  y  for  gasses  is  the  pressure  de¬ 
pendence.  y  is  proportional  to  the  number  density  of  the  molecules  in  the  gas.  This  can 
be  seen  from  the  dispersion  equation, 


«2(o>)  =  1  + 


tjne 


(5.7) 


where  n  is  the  index  of  refraction,  aj  is  the  frequency  of  the  light,  A'  is  the  number  of 
molecules  per  unit  volume  in  the  gas,  q,  is  the  electron  charge,  m,  is  the  electron  mass, 
f.  is  the  number  of  oscillators  in  the  jth  molecule,  and  oj0,  is  the  resonance  frequency  of 
the  jth  molecule.  [Ref.  10:  p.  60]  Defining  the  refractivity, 

X  =  n  -  1,  (5.8) 

the  first  order  Taylor  expansion  of  the  left  hand  side  of  (5.7)  is 

1  +  2y.  (5.9) 


This  approximation  is  valid  for  C02  since  y<l  .  Substituting  (5.9)  into  (5.7)  yields 


Xoc.Y,  (5.10) 

the  number  density  of  molecules  in  the  gas.  For  an  ideal  gas,  the  relationship  between 
number  density  and  the  quantities,  pressure  and  temperature,  is  specified  by 


A'  = 


P 

RT  ' 


(5.11) 


where  A'  is  the  number  density,  R  is  the  ideal  gas  constant,  and  T  is  the  temperature  of 
the  gas.  [Ref.  11:  p.  470]  Thus  the  number  density 


This  leads  to 


(5.12) 
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Both  of  the  modified  codes  utilize  this  relationship  to  calculate  x.  for  the  gasses.  A 
table  of  values  for  indicies  of  refraction  at  standard  atmospheric  pressure  and  at  293° 
kelvin  are  used  as  a  basis  for  computation.  Values  for  pressure  and  temperature  are  read 
in  from  the  users  input  file  and  used  to  scale  x •  This  eliminates  the  need  for  the  user  to 
enter  a  new  table  for  index  of  refraction  for  the  gasses  ever}'  time  he  desires  to  investigate 
the  behavior  at  a  different  pressure  or  temperature. 

D.  SENSITIVITY  ANALYSIS 

To  test  the  sensitivity  of  the  codes  to  variation  in  some  of  the  key  parameters,  runs 
were  made  using  different  indicies  of  refraction  and  different  beam  energies.  A  run  was 
made  using  carbon  dioxide  at  standard  atmospheric  pressure  as  opposed  to  2  atmos¬ 
pheres,  as  in  the  previous  runs.  At  1  atmosphere  of  C02,  the  indicies  of  refraction  for 
C02  are  below  the  cutoff  for  Cerenkov  production  with  the  electron  beam  energy  at  16.7 
MeV.  (The  Cerenkov  threshold  at  1800  angstroms  and  1  atmosphere  is  14.9  MeV.)  This 
means  that  ideally  there  is  no  production  of  Cerenkov  radiation  in  the  carbon  dioxide 
(zone  5).  (Other  sources  such  as  transition  radiation  have  not  been  addressed  in  this 
work.)  An  increase  in  beam  energy  again  causes  np  >  1  and  produces  Cerenkov  radi¬ 
ation  in  the  C02  (zone  5).  Note  that  in  Figure  22  on  page  41,  there  is  no  Cerenkov 
production  in  zone  5  for  the  16.7  MeV  beam  while  there  is  for  the  17  Mev  beam.  Both 
the  16.7  MeV  and  17  MeV  beams  produced  about  the  same  amount  of  Cerenkov  radi¬ 
ation  as  the}  passed  through  the  glass  mirror  (zone  6). 

Since  the  indicies  of  refraction  of  gasses  are  dependent  upon  the  pressure  of  the  gas. 
the  modified  codes  were  designed  to  .allow  the  user  to  enter  the  pressures  of  any  air  or 
carbon  dioxide  in  the  geometry.  The  codes  then  calculate  the  indicies  of  refraction  for 
these  gasses.  The  pressure  and  temperature  of  the  C02  and  of  the  air  are  a  part  of  the 
regular  input  file.  (Appendix  E). 

These  results  show  that  the  output  is  very  sensitive  to  the  indicies  of  refraction  used 
in  the  code  patch  and  to  the  energy  of  the  source  particle  beam.  This  is  important  since 
it  allows  the  user  to  vary  parameters  in  the  problem  and  test  the  results.  For  example, 
a  user  could  experiment  with  different  Cerenkov  detector  designs  using  different  materi¬ 
als  and  thicknesses  for  the  converter  foil,  different  pressures  of  C02,  or  try  glass  or  quartz 
to  produce  the  Cerenkov  radiation.  He  may  also  be  able  to  select  certain  regions  of  the 
spectrum  to  emphasize  based  upon  how  the  Cerenkov  light  is  refracted  through  his  op¬ 
tics. 
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Variations  in  Pressure  and  Energy:  This  figure  shows  the  Cerenkov 
production  in  C02  at  1  atmosphere  by  a  16.7  MeV  electron  beam  and 
by  a  17  MeV  electron  beam. 
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E.  FIGURE  OF  MERIT 


The  Figure  of  merit  of  a  Monte  Carlo  computer  model  is  defined  as 


where  a  represents  the  1  a  relative  error  and  /  is  the  computer  run  time  to  achieve  the 
same  error.  The  utility  of  the  figure  of  merit  is  that  it  allows  the  user  to  compare  the 
relative  efficiencies  of  different  methods  used  on  a  computer  (the  larger  the  better).  The 
figure  of  merit  was  calculated  on  the  IBM  3033  computer.  In  the  detector  zone  the  fig¬ 
ure  of  merit  for  the  MT  method  was  2.71  x  10-3.  This  was  uniform  across  the  entire 
range  of  wavelengths.  The  figure  of  merit  for  the  DS  method  was  3.37  x  10-3  for  wave¬ 
lengths  shorter  than  5000  angstroms  and  was  2.14  x  10-3  for  wavelengths  of  5000 
angstroms  and  longer.  Thus,  for  wavelengths  shorter  than  5000  angstroms,  the  DS  fig¬ 
ure  of  merit  was  24  percent  higher  than  the  MT  figure  of  merit.  However,  for  wave¬ 
lengths  longer  than  5000  angstroms,  the  MT  figure  of  merit  was  27  percent  higher  than 
the  DS  figure  of  merit.  This  problem  with  the  DS  method  can  be  removed  by  forcing 
DS  to  sample  the  longer  wavelengths  more  often  provided  that  the  results  of  this  biasing 
are  removed  by  a  corresponding  weight  adjustment  (the  weight  assigned  to  the  longer 
wavelengths  would  have  to  be  reduced  and  the  weight  assigned  to  shorter  wavelengths 
increased  to  insure  proper  weight  balance). 

While  these  figures  of  merit  tend  to  support  the  idea  that  the  best  code  for  the 
computer  time  available  is  the  distribution  sampling  method,  the  user  should  be  cau¬ 
tioned  that  this  is  only  true  for  wavelengths  shorter  than  5000  angstroms.  The  sampling 
error  for  the  longer  wavelengths  causes  this  method  to  lose  value  to  anyone  interested 
in  the  entire  Cerenkov  spectrum  from  1800  to  10000  angstroms.  However,  since  this 
code  performs  well  in  the  region  where  the  most  Cerenkov  production  takes  place,  it  is 
a  good  choice  for  anyone  interested  in  holding  run  time  down  to  a  minimum  while  still 
getting  wavelength  dependent  information. 

The  Cl  code  patch  is  the  fastest  of  all  the  methods  but  predicts  too  little  Cerenkov 
production  at  the  shorter  wavelengths  and  too  much  at  the  longer  wavelengths.  It  is  still 
valid  for  the  user  who  is  not  interested  in  detailed  wavelength  dependent  information 
about  the  Cerenkov  spectrum.  The  MT  routine  produces  the  most  accurate  and  detailed 
results  of  the  three  code  patches;  however,  it  requires  substantial  computer  time. 
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VI.  CONCLUSIONS 


The  constant  index  (Cl)  Cerenkov  routine  provides  a  valuable  tool  in  the  design  of 
Cerenkov  detectors  and  for  the  study  of  Cerenkov  radiation.  However,  the  use  of  a 
constant  index  of  refraction  in  the  calculation  causes  the  routine  to  under-estimate  the 
Cerenkov  production  at  short  wavelengths  and  to  over-estimate  production  at  longer 
wavelengths.  For  carbon  dioxide  gas,  this  produces  an  error  in  the  Cl  production 
weights  of  -44  percent  at  1800  angstroms  and  +5  percent  at  10000  angstroms.  These 
figures  are  based  on  trial  runs  using  n  =  1.000894  as  the  constant  index  of  refraction  in 
the  Cl  method.  Even  if  a  user  selects  a  constant  index  of  refraction  so  that  the  total 
Cerenkov  weight  produced  is  the  same  for  Cl  as  for  the  DS  and  MT  patches,  the  Cl 
method  will  under-estimate  the  Cerenkov  production  in  the  iSuO-angstrom  bin  by  about 
23  percent  while  over-estimating  the  production  in  the  10000-angstrom  bin  by  about  8 
percent.  This  discrepancy  is  caused  by  the  sensitivity  of  Cerenkov  production  to  the 
refractivity  (n  -  1)  of  the  material  in  which  production  occurs. 

The  Cl  approach  has  been  modified  in  the  multiple  trace  (MT)  method.  MT  uses  a 
brute  force  approach  to  include  wavelength  dependence  and  provides  the  most  detailed 
treatment  of  the  three  code  patches.  It  allows  the  user  to  study  spectral  details  not 
available  in  the  Cl  method  but  at  a  significant  cost  in  computer  time. 

A  compromise.  the  distribution  sampling  (DS)  method,  has  also  been  provided. 
This  approach  has  a  figure  of  merit  that  is  25  percent  higher  than  the  MT  method  for 
wavelengths  shorter  than  5000  angstroms.  However,  the  figure  of  merit  drops  below 
that  of  MT  for  wavelengths  longer  than  5000  angstroms  due  to  the  small  number  of 
times  these  wavelengths  are  sampled.  Since  80  percent  of  the  Cerenkov  production  oc¬ 
curs  below  5000  angstroms,  this  patch  can  be  used  effectively  by  arnone  that  does  not 
need  information  on  the  entire  spectrum.  The  problem  with  DS  at  the  longer  wave¬ 
lengths  can  be  removed  by  biasing  the  sampling  so  that  all  wavelengths  are  sampled 
uniformly.  However,  care  must  be  taken  to  adjust  the  weights  so  that  the  weight  as¬ 
signed  to  the  longer  wavelength  bins  is  reduced  while  the  weight  assigned  to  the  shorter 
wavelength  bins  is  increased  in  order  to  maintain  proper  weight  balance. 

The  MT  and  DS  code  patches  yield  additional  detail  on  the  Cerenkov  spectrum 
generated  by  energetic  electrons.  Together  with  the  Cl  routine,  the  MT  and  DS  mod¬ 
ifications  to  ITS  provide  the  user  with  an  excellent  set  of  tools  for  use  in  the  study  of 
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Cerenkov  radiation  and  in  the  design  of  related  diagnostic  equipment.  There  remains  the 
requirement  to  have  a  thorough  verification  against  experimental  results. 
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APPENDIX  A.  ITS  (ACCEPT) 


ACCEPT  is  the  ITS  code  for  general,  combinatorial  geometries.  This  addition  of 
Cerenkov  radiation  to  the  ITS  code  works  only  with  the  ACCEPT  code. 

To  run  ITS  (ACCEPT)  the  user  is  required  to  do  the  following: 

1.  Create  a  compiled  version  of  the  code  including  the  subroutine  changes  to  the  basic 
code.  This  is  done  by  running  the  UPEML  program  pros'ided  with  the  ITS  pack- 
aee.  The  input  files  to  this  L’PEML  run  are  provided  as  Appendix  B  and  Appendix 
C. 

2.  Create  a  cross  section  tape  by  running  the  program  XGEX  which  is  a  part  of  the 
ITS  package.  An  example  input  file  is  provided  in  Appendix  D. 

3.  Execute  the  modified  ACCEPT  program.  An  example  problem  input  file  is  also 
provided  in  Appendix  E  along  with  a  description  of  the  combinatorial  geometry 
scheme. 
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APPENDIX  B.  MULTIPLE  TRACE  SUBROUTINE  L’  liNG 

*IDENT  DEFINE 
*DEFINE  IBM 
^DEFINE  DOUBLE 
*DEFINE  ACCEPT 
* I DENT, PHILLIPS 
*IDENT,PERMT 
*D,PARAMS.  17 

PARAMETER  (  KPTMAX=10000,  INSTAT=100, 

*1 ,HIST.  94 

COM  =  CZERO 
*I,HIST.  223 

IF( ( IB  .EQ.  1)  .AND.  flHIST  . LE.  50))  THEN 
IF( IHIST  .EQ.  1)  WRITE( 6 , 801 ) 

C  SEND  ADITIONAL  OUTPUT  TO  OUTPUT  UNIT  6 

801  FORMAT ( IX , 1HI , 2X , 2HLB , 7X , 1HX , 14X , 1HY , 14X , 1HZ , 

1  13X, 1HU, 13X, 1HV, 12X, 1HV, 12X, 3HCOM, 10X, 7HENERGY  /) 

WRITE( 6 , 802)  IHIST,LB,XB,WT,COM,T 
802  FORMAT( 213 , 8( 2X , 1PE12. 5) ) 

C 

ENDIF 

C 

*IDENT,CKVTR 
*1 ,CALC.  78 

6  ,AIRP,C02 ,TMPKEL 
*I,OUT.  30 

C  CERENKOV  ADDITIONS  TO  COMMON  BLOCK 
COMMON  /OUT/ 

1  CHER(  5 , 100) , CHRNUM( 100 , 100) , 

2  CHRWAT( 100 , 100) ,CHRENG( 100 , 100) ,CHRABS( 100,100) , 

3  CHROUT( 100 , 100) ,CHRENT( 100 , 100) ,TRNUM( 100 , 100) , 

4  TRWAT( 100, 100) ,TRENG( 100, 100) ,TRABS( 100 , 100) , 

5  TROUT ( 100,100) ,TRENT( 100 , 100) ,NRVIR ,NCVIR ,NAVIR , 

6  CVIRB( 11) .AVIRBC 11) ,RVIRB( 11) ,VARIND( 100,100) ,WVLNTH(100) , 

7  S  CO VOL( 1 0 , 1 0 , 1 0 ) , HV I R , I TR AD , NRB  OTH , I CKV , I V I RS , KALTOT , 

8  LBTRAD , NBOTRD , JSRTRD , NCALL , NMB INS , WKILL 
*1 , INPUT. 54 

DIMENSION  WVLNT( 100) ,VARIN( 100,100) 


C*****THIS  SECTION  ASSIGNS  VALUES  TO  THE  TABLE  OF  INDEX  OF  ***** 
C*****REFRACTION  AND  WAVELENGTH.  THE  TABLES  CREATED  ARE  ***** 

C*****STRUCTURED  TO  COVER  WAVELENGTHS  IN  200  ANGSTROM  BINS  ***** 

C*****FROM  1800  THROUGH  5000  ANGSTROMS,  AND  IN  500  ANGSTROM  ***** 
C*****BINS  FROM  5500  THROUGH  10000  ANGSTROMS.  ***** 

C*****THE  VARIND  TABLE  CONTAINS  INDICES  OF  REFRACTION  FOR  ***** 

C*****EACH  MATERIAL.  ***** 

C*****  ***** 

C*****VARIND( I , J)  --  INDEX  OF  REFRACTION  FOR  MATERIAL***** 

C*****  MATERIAL  I  ASSOCIATED  WITH  ***** 

C*****  WAVELENGTH  BIN  J  ***** 

(^-.VVrjWrfr  ^WnY'jV'sV 

C*****WVLNTH(J)  --  WAVELENGTH  VALUE  FOR  BIN  J  ***** 
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Q-jV'jVtV’jV’jV  ickick'k 

c 

C  WAVELENGTHS 

DATA  (WVLNT( I), 1=1,27)/  1800,2000, 

1  2200,2400,2600,2800,3000, 

2  3200,3400,3600,3800,4000, 

3  4200,4400,4600,4800,5000, 

4  5500,6000,6500,7000,7500, 

5  8000,8500,9000,9500, 1G00G  / 

C  VACUUM 

DATA  ( VARIN( 1,I),I=1,27)/27*1. DO/ 

C  AIR 

DATA  (VARIN(2,I) ,1=1,27)/ 

1  1. 00034D0 , 1. 0003100,1. 00030D0.1. 00030D0,1. 00029D0, 

2  1.  0002900,1. 00029D0 , 1. 0002900,1. 00028DO,1. 00028D0, 

3  1. 00028D0 , 1. 0002800,1. 00028D0,1. 00028D0,1. 00028D0, 

4  1. 00028D0 , 1.  00028DO , 1. 00028D0,1. 00028D0.1. 00028D0, 

5  1. 00028D0 , 1. 0002700,1. 00027D0,1. 00027D0.1. 00027D0, 

6  1. 00027D0 , 1. 00027DO  / 

C  CARBON  DIOXIDE  (STANDARD  ATMOSPHERE) 

DATA  (VARIN(3,I) ,1=1,27)/ 

1  1. 00055D0,1.  00052D0,1.  00051D0,1. 00050D0.1. 00049D0, 

2  1.  00048D0 , 1.  00048D0 , 1. 00047D0,1. 00047D0.1. 00046D0, 

3  1. 00046D0 , 1. 00046D0 , 1. 00046D0,!. 00046D0.1. 00045D0, 

4  1.  00045D0 , 1.  00045D0 , 1. 00045D0.1. 00045D0,1. 00045D0, 

5'  1  00045D0 , 1. 00045D0 , 1. 00045D0,1. 00044D0.1. 00044D0, 

6  1. 00044D0 , 1. 00044D0  / 

C  GLASS  TYPE  458/678(MIL-G- 174)  (FUSED  SILICA  DENSITY  2.202  G/CM**3) 
DATA  (VARIN(4, I) ,1=1,27)/ 

1  1. 58529D0 , 1. 5505100,1. 52845D0.1. 51333D0.1. 50352D0, 

2  1.  49565D0 , 1. 48779D0 , 1. 48274D0 , 1. 47884D0,1. 47554D0, 

3  1. 47283D0 , ] . 47012D0 , 1. 46830D0.1. 46648D0,1. 46492D0, 

4  1.  4356800, .  46233D0 , 1.  4599  IDO,  1.  45804D0,1,  45664D0, 

5  1. 4552*D0, 1. 4542400,1. 45332D0.1. 45250D0, 1. 45175D0, 

6  1. 45107D0, 1. 45042D0  / 

C  CRYSTAL  QUARTZ 

DATA  (VARIN(5, I) ,1=1,27)/ 

1  1. 6600600,1. 6508700,1. 62626D0.1. 61011D0,1. 60158D0, 

2  1. 59306D0 , 1. 58453D0, 1. 5760000,1. 56747D0, 1. 56413D0, 

3  1. 56080D0 , 1. 5577900,1.  5555400,1. 55349D0.1.  55194D0, 

4  1. 55039D0 , 1 . 54884D0.1. 54616D0.1. 5439300,1. 54247D0, 

5  1. 54 10  IDO, 1. 5  3955D0, 1. 5383900,1. 5374400,1. 5366300, 

6  1. 53581D0 , 1. 5350200  / 

C 

*1 .INPUT. 436 

C  CERENKOV  ADDITIONS  TO  DEFAULT  INPUT  PARAMETERS 
NMBINS=27 
NCALL=0 
ICKV=0 
IVIRS=0 
ITRAD=0 

C  END  OF  CERENKOV  ADDITIONS 
*1 .INPUT. 575 

C ********************** 
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C  CERENKOV 

C********************** 

ELSE  IF  (KOP( 'CERENKOV' )  .  GE.  0)  THEN 
IMATCH  =  1 
ICKV  =  1 

C  SET  PRESSURE  FOR  AIR  IN  STANDARD  ATMOSPHERES 
AIRP  =  PARM(l) 

C  SET  PRESSURE  FOR  C02  IN  STANDARD  ATMOSPHERES 
C02P  =  PARM( 2) 

C  SET  THE  TEMPERATURE  IN  KELVIN 
TMPKEL  =  PARM( 3) 

C  SCALE  TEMPERATURE  TO  293  DEGREES  KELVIN 
TMPKEL  =  TMPKEL  /  293. DO 
DO  93774  I 1=1 ,NMBINS 

VARIND( 1 , I1)=VARIN( 1,11) 

C  INDICES  OF  REFRACTION  -  1  CHANGE  PROPORTIONAL  TO  PRESSURE  CHANGE 
VARIND( 2,11 )=( VARIN( 2 , II) -1.  DO)*AIRP/TMPKEL  +  l.DO 
VARIND(3, 11)=(VARIN(3, II) -1.  D0;*CO2P/TMPKEL  +  l.DO 
93774  CONTINUE 

DO  93884  Jl=4 , 20 

DO  93874  1 1=1 ,NMBINS 

VARIND( J1 , I 1)=VARIN( J1 ,11) 

93874  CONTINUE 

93884  CONTINUE 

DO  93894  I 1=1 ,NMBINS 

WVLNTHC I 1 )=WVLNT( I 1 ) 

93894  CONTINUE 

KALTOT  =  PARM(4) 

DO  3034  JJJ=1 ,NZON 

CALL  OPREADC 1 , INUNIT, IECHO, IVAL,EOFLAG) 

C 

C  CHECK  INPUT  FILE  FOR  COMPLETE  OPTICAL  DATA 
C 

IF(EOFLAG)  THEN 

WRITE( * ,*)  'PREMATURE  END  OF  OPTICAL  DATA' 

STOP 
END  IF 

IF  (KOP( ’+’ )  . GE.  0)  THEN 
C 

C  ASSIGN  VALUES  TO  CERENKOV  INPUT  MATERIAL  NUMBER  FOR  ZONE  JJJ 

CHER( 1 , JJJ)  =  PARM(l) 

C  EXTINCTION  COEFFICIENT  (1/CM)  --  ATTENUATION  COEFFICIENT 

CHER( 2, JJJ)  =  PARM( 2) 

C  FRACTION  OF  LIGHT  TO  BE  MIRROR  REFLECTED  FROM  ZONE  WALLS 

CHER( 3 , JJJ)  =  PARM( 3) 

C  FRACTION  OF  LIGHT  TO  BE  DIFFUSE  REFLECTED  FROM  ZONE  WALLS 

CHER(4 , JJJ)  =  PARM( 4) 

C  FRACTION  OF  LIGHT  TO  BE  ACCEPTED  FOR  ENTRANCE  INTO  ZONE 

CHER(5 , JJJ)  =  PARM( 5) 

C 

ELSE 

WRITE(*,8903) 

8903  FORMAT( IX, ’ ILLEGAL  CKV  OR  VIRS  ZONE  DATA’) 

GO  TO  9999 
END  IF 

3034  CONTINUE 
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0********************* 

C  VIRTUAL-SOURCE 

ELSE  IF  (KOP( 'VIRTUAL-SOURCE' )  . GE.  0)  THEN 
I MATCH  =  1 
IVIRS  =  1 
HVIR  =  PARM(l) 

NRVIR  =  PARM( 2) 

NCVIR  =  PARM( 3) 

NAVIR  =  PARM(4) 

READ( INUNIT, * ,END=905 1 ,ERR=9052) 

1  (RVIRB(IR) ,IR=1 .NRVIR)  ,RVIRB(NRVIR  +  I) 

READ( INUNIT,* ,END=9051 ,ERR=9052) 

1  (CVIRB( IC) , IC=1 , NCVIR) ,CVIRB( NCVIR  +  1) 

READ( INUNIT, *,END=905 1 ,ERR=9052)( AVIRB( IA) , 

1  IA=1, NAVIR)  ,AVIRB(NAVIR  +  1) 

q ********************** 

C  TRANSITION  RADIATION 

QVoWr  ‘V  iV  •>>  Vr  -5 V  *  iV  Vc  *  *  *  *  *  Vr  Vc  Vr  Vc *  Vr 

ELSE  IF(K0P( 'TRANSITION-RADIATION' )  .  GE.  0)  THEN 
I MATCH  =  1 
ITRAD  =  1 
LBTRAD  =  PARM(l) 

NBOTRD  =  PARM( 2) 

JSRTRD  =  PARM( 3) 

*1, INPUT. 1178 

IF( ICKV  .EQ.  1)  THEN 

C  CERENKOV  OPTICAL  AND  VIRTUAL  SOURCE  DATA . 

WRITE( 6 , 9005) 

9005  FORMATC 6H  Z0NE,6X,6H  MTRL  , 6X , 6HABS0RB , 6X , 6HMIRROR , 5X, 7HDIFFUSE 
1  , 6X , 6HACCEPT , 4X , 8HTRANSMIT) 

DO  90051  1=1 , LZMAX 

C  FRACTION  OF  LIGHT  TO  BE  TRANSMITTED  OUT  THROUGH  ZONE  BOUNDARY 
FST  =  1  -  CHER( 3,1)  -  CHER(4,I) 

WRITE( 6 , 90050) I , ( CHER( J , II , J=1 ,5) ,FST 

90050  FORMATC 16 , 1P6E12.  4) 

90051  CONTINUE 
C 

END  IF 
C 

IF( ITRAD  .EQ.  1)  THEN 

WRITE( 6 ,90311)  LBTRAD .NBOTRD , JSRTRD 
90311  FORMATC/'  TRANSITION  RADIATION  LOCATION,  ZONE  BODY  SURFACE' 

1  ,316) 

C 

ENDIF 

C 

IF  (IVIRS  .EQ.  1)  THEN 

WRITE( 6,90300)  HVIR , NRVIR , NCVIR , NAVIR 
90300  FORMATC /17H  VIRTUAL  SOURCE  Z,  E12.4 

1  /'  VIRTUAL  SOURCE  TALLY  BINS,  RADIAL,  POLAR,  AZIMUTH,' 

2  ,  316) 

WRITE(  6 ,9050)  (RVIRB( IR) , IR=1 , NRVIR) ,RVIRB( NRVIR  +  1) 

WRITE ( 6,9055)(CVIRB(IC) ,IC=1 , NCVIR) ,CVIRB( NCVIR  +  1) 

WRITE( 6 ,9060)( AVIRB( IA) , IA=1 , NAVIR) ,AVIRB( NAVIR  +  1) 

9050  FORMATC 8H  RADIAL  .10E12.4) 
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9055  FORMAT(8H  POLAR  .10E12.4) 

9060  F0RMAT(8H  AZIMUTH  ,10E12. 4) 

C 

END  IF 
C 

*1 ,PREP.  659 

IFCICKV  .EQ.  1)  THEN 
WKILL  =  CZERO 
DO  241  LLL=1 , 100 
DO  2401  JJJ=1,100 

CHRNUMC LLL , JJJ )  =  CZERO 
CHRWATC LLL, JJJ)  =  CZERO 
CHRENG(LLL, JJJ)  =  CZERO 
CHRABS(LLL, JJJ)  =  CZERO 
CHROUTC LLL , JJJ)  =  CZERO 
CHRENTC LLL , JJJ )  =  CZERO 
2401  CONTINUE 

241  CONTINUE 
C 

ENDIF 

C 

IF( IVIRS  .EQ.  1)  THEN 
DO  244  IR=1 ,NRVIR 
DO  243  IC=1 ,NCVIR 
DO  242  IA=1 ,NAVIR 

SCOVOL( IA , IC , IR)  =  CZERO 

242  CONTINUE 

243  CONTINUE 

244  CONTINUE 
C 

ENDIF 

C 

*1, OUTPUT.  47 

DIMENSION 

1  CHRN( 100) ,CHRW( 100) ,CHRE( 100) ,CHRA( 100) ,CHRO( 100) ,CHRI( 100) 
INTEGER  ISIGA( 100) , ISIGB( 100) , ISIGCC 100) , ISIGD( 100) , ISIGE( 100) , 

1  ISIGF( 100) , IWVL( 100) , ISIG1( 100)  ,ISIG2(100)  ,ISIG3(100)  , 

2  ISIG4( 100)  , ISIG5( 100)  ,ISIG6(100) 

*1 .OUTPUT.  561 

IFCICKV  .EQ.  1)  THEN 
WTENT  =  CZERO 
WTABS  =  CZERO 
WTOT  =  CZERO 
WTCON  =  CZERO 
DO  66  1=1 , LZMAX 

DO  56110  J=1 ,NMBINS- 1 

CHRNUMC I , J)=(CHRNUM( I , J)/(CIMAX) ) 

CHRWATC I , J)=CHRWAT( I , J)/CIMAX 
CHRENGC I , J)=CHRENG( I , J)/CIMAX 
CHRABSC I , J)=CHRABS( I , J)/CIMAX 
CHROUTC I , J ) =CHROUT( I , J ) / C IMAX 
CHRENTC I , J)=CHRENT( I , J) /CIMAX 
56110  CONTINUE 

66  CONTINUE 

DO  10225  1=1,100 
CHRN( I )=CZERO 
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CHRW( I )=CZERO 
CHRE( I)=CZERO 
CHRA( I )=CZERO 
CHRO( I)=C2ER0 
CHRI( I )=CZERO 
10225  CONTINUE 

C  COMPUTE  THE  TOTALS  FOR  EACH  ZONE 
DO  11230  1=1 , LZMAX 

DO  11225  J=1 ,NMBINS-1 

CHRN ( I )  =CHRN(  I  )■ +CHRNUM ( I ,  J ) 

CHRW( I )=CHRW( I )+CHRWAT( I , J) 

CHRE( I )=CHRE( I )+CHRENG( I , J) 

CHRA( I )=CHRA( I )+CHRABS( I , J) 

CHRO( I)=CHRO( I )+CHR0UT( I , J) 

CHR I ( I ) =CHR I ( I ) +CHRENT( I , J ) 

11225  CONTINUE 

11230  CONTINUE 

C 

C  WEIGHT  CONSERVATION  FRACTION  =  (TOTAL  WEIGHT  ABSORBED  +  TOTAL 
C  WEIGHT  KILLED  +  WEIGHT  THAT  ESCAPES)  /  WEIGHT  PRODUCED. 

WKILL  =  WKILL/CIMAX 
WTENT  =  CHRI(LZMAX) 

DO  13623  1=1 , LZMAX 

WTABS  =  WTABS  +  CHRA(I) 

WTOT  =  WTOT  +  CHRW(I) 

13623  CONTINUE 

WTCON  =  (WKILL  +  WTABS  +  WTENT) /WTOT 
END  IF 
C 

IF  ( IVIRS  .EQ.  1)  THEN 
DO  96  IR=1 ,NRVIR 
DO  86  IC=1 ,NCVIR 
DO  76  IA=1,NAVIR 

SCOVOL( IA , IC , IR)=SCOVOL( I A , IC , IR) /CIMAX 
76  CONTINUE 

86  CONTINUE 

96  CONTINUE 

C 

ENDIF 

C 

*1, OUTPUT.  1375 

IF( ICKV  .EQ.  1)  THEN 

C********Vr*********cj{^;Q£s  ADDED  here********************* 
NPUT=NMB INS - 1 
TEMP( 1)=WTC0N 

C  CALCULATE  AVERAGES  AND  STATISTICS  FOR  TOTALS 
CALL  STATS (TEMP, ISIG, 1 ,KPUT) 

WTCON=TEMP( 1) 

CALL  STATS( CHRN, ISIG1, LZMAX ,KPUT) 

CALL  STATS (CHRW, IS IG2, LZMAX, KPUT) 

CALL  STATS(CHRE,ISIG3, LZMAX, KPUT) 

CALL  STATS ( CHRA , IS IGA , LZMAX , KPUT) 

CALL  STATS ( CHRO , I S I G5 , LZMAX , KPUT ) 

CALL  STATS ( CHR I , I S I G6 , LZMAX , KPUT ) 

DO  11220  1=1, LZMAX 
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C  COMPUTE 


11110 


11115 


11120 


11125 


11130 


11135 


11140 


11145 


11150 


11155 


11160 


11165 


11170 


WRITE(NPRT,9000)  I 
FORMAT( / , ' ZONE  =  '  ,13,/, 

1 WVLNTH 1 ,T14, ’NUMBER’ ,T22,’STAT’ ,T33 ,’ WEIGHT’ ,T41,’STAT’ , 
T52, ’ENERGY' ,T60, ’STAT’ ,T7 1 ,' ABSORB ’ ,T79, ’STAT’ ,T90, 

’  LEAVE’ ,T97,’STAT’ ,T108,’  ENTER’ ,T1 15 ,’ STAT* ,/ , 

’ . ’  ,T14 ,  ’ . '  ,T22,  ’ - '  ,T33 , ' . ’  ,T41,  ’ - ’  , 

T52 ,  ’ . ’  ,T60 ,  ’ - ’  ,T7 1 ,  ’ . ’  ,T79,’ - ’  ,T90, 

’ . ’  ,T97  ,  ’ - ’  ,T108,  ’ . ’  ,T115  ,  ’ - ’  ,/) 

THE  AVERAGES  AND  THE  PERCENT  ERROR 
DO  11110  J=1 .NMBINS-l 
TEMP( J)=CHRNUM( I , J) 

CONTINUE 

CALL  STATS ( TEMP , I S I GA , NPUT , KPUT ) 

DO  11115  J=1 ,NMBINS- 1 
CHRNUM( I , J)=TEMP( J) 

CONTINUE 

DO  11120  J=1,NMBINS-1 
TEMP( J)=CHRWAT(  I , J) 

CONTINUE 

CALL  STATS(TEMP,ISIGB, NPUT, KPUT) 

DO  11125  J=1 , NMBINS-l 
CHRWAT( I , J)=TEMP( J) 

CONTINUE 

DO  11130  J=1 .NMBINS- 1 
TEMP( J)=CHRENG( I , J) 

CONTINUE 

CALL  STATS( TEMP, I SIGC, NPUT, KPUT) 

DO  11135  J=1 ,NMBINS-1 
CHRENG( I , J)=TEMP( J) 

CONTINUE 

DO  11140  J=1 ,NMBINS-1 
TEMP( J)=CHRABS( I , J) 

CONTINUE 

CALL  STATS(TEMP,ISIGD, NPUT, KPUT) 

DO  11145  J=1 .NMBINS- 1 
CHRABS( I , J)=TEMP( J) 

CONTINUE 

DO  11150  J=1 .NMBINS- 1 
TEMPC J)=CHROUT(  I , J) 

CONTINUE 

CALL  STATS(TEMP,ISIGE, NPUT, KPUT) 

DO  11155  J=1 .NMBINS-l 
CHROUT( I , J)=TEMP( J) 

CONTINUE 

DO  11160  J=l, NMBINS-l 
TEMP( J)=CHRENT( I , J) 

CONTINUE 

CALL  STATSC TEMP, IS IGF, NPUT, KPUT) 

DO  11165  J=l, NMBINS-l 
CHR£NT( I , J)=TEMP(  J) 

CONTINUE 
DO  11170  J= 1,17 

IWVL( J)=1600  +  J*200 
CONTINUE 

DO  11180  J=18 , 27 
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11181 


11190 

11220 

90055 


IWVL(J)=5000  +  ( J- 17 )*500 

11180  CONTINUE 

DO  11190  J=1 .NMBINS-l 

WRITECNPRT, 11181)  ItfVL( J) ,CHRNUM( I ,J) ,ISIGA( J) , 

1  CHRWAT(I.J), 

2  ISIGBC J) ,CHRENG( I ,  J) ,  ISIGC( J) ,CHRABS( I , J) ,ISIGD( J) , 

3  CHROUT(I.J) ,ISIGE( J) ,CHRENT(I,J) ,ISIGF( J) 

11181  FORMAT ( 15 ,T8 , 1PE12.  4 ,T22 , 13 ,T27 , 1PE12.  4  ,T41 , 13 ,T46 , 

1  1PE12.4, 

2  T60 , 13 ,T65 , 1PE12. 4 ,T79 , 13 ,T84, 1PE12. 4 ,T97 , 13 ,T102 , 

3  1PE12. 4,T115 , 13) 

11190  CONTINUE 

11220  CONTINUE 

WRITE(NPRT, 90055) 

90055  FORMAT! / . ' TOTALS  FOR  EACH  ZONE',/, 

1  ’ ZONE* ,T14, 'NUMBER' ,T22  'STAT' ,T33 , 'WEIGHT' ,T4l ,'STAT’ , 

2  T52, 'ENERGY' ,T60, ' STAT' ,T71 , 'ABSORB' ,T79 , ' STAT' ,T90, 

3  ’  LEAVE' ,T97, 'STAT' ,T108, :  ENTER’ ,T115 ,' STAT' ,/ , 

1  ' . '  ,T14 , ' . ’  ,T22  , ' - '  ,T33,' . ’  ,T41,' - 

2  T52  , ' . '  ,T60, ' - '  ,T71,  ' . '  ,T79,  ' - ’  ,T90, 

3  ' . ’  ,T97  , ' - '  ,T108,  ' . ’  ,T115  ,  ' - '  ,/) 

DO  17190  J=1 , LZMAX 

WRITECNPRT, 17181)  J,CHRN( J) , ISIG1C J) , 

1  CHRW( J) ,ISIG2(J)  , 

2  CHRE( J) , ISIG3( J)  ,CHRA( J) ,ISIG4( J) , 

3  CHRO(J) ,ISIG5( J) ,  CHRI(J),ISIG6(J) 

17181  F0RMAT(I5,T8,1PE12.  4 ,T22 , 13 ,T27 , 1PE12.  4 ,T41 , 13 ,T46 , 

1  1PE12.4, 

2  T60 , 13 ,T65 , 1PE12.  4 ,T79 , 13 ,T84 , 1PE12. 4 ,T97 , 13 ,T102 , 

3  1PE12. 4.T115 , 13) 

17190  CONTINUE 

WRITECNPRT, 11246)  WTCON 

11246  FORMAT( / , '  WEIGHT  CONSERVATION  FRACTION  IS  ' ,E12. 4) 

END  IF 

C**********************£ND  OF  CHANGED  OUTPUT*********************** 
90100  FORMAT! I3,3X,6(E12.  4,13)) 

IF( IVIRS  .EQ.  1)  THEN 
WRITECNPRT, 90200) 

90200  FORMAT! /' CERENKOV  VIRTUAL  SOURCE,  PHOTONS / SQCM/ STR ' ) 

DO  6030  IVR=1 ,NRVIR 

WRITECNPRT, 90300)  RVIRBCIVR) ,RVIRB( IVR+1) 

90300  FORMAT! '  R-MIN' ,E12. 4/ '  R-MAX' ,E12. 4,24X, 

1  '  UPPER  AND  LOWER  POLAR  COSINES  ’) 

DO  6020  LVC=1 ,NCVIR ,5 
MVC=MINO( LVC+4.NCVIR) 

WRITECNPRT , 9040 ) ( CVIRB ( I VC ) , I VC=LVC , MVC ) 

WRITECNPRT, 9050) ( CVIRBC IVC+1 ) , IVC=LVC ,MVC) 

9040  FORMAT! 23X, 'MINIMUM' ,5X, ' MAXIMUM* ,6E15.  4) 

9050  FORMAT! 2 3X, 'AZIMUTH' ,5X, 'AZIMUTH' ,6E15.  4) 

DO  6010  IVA=1 ,NAVIR 
L=0 

DO  6000  IVC=LVC,MVC 
L=L+1 

TEMP! L)=SCOVOL( IVA , IVC , IVR) 

6000  CONTINUE 

CALL  STATS(TEMP , ISIG ,L,KPUT) 


17181 


17190 


11246 


90200 


90300 
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WRITE( NPRT , 9060 ) AVIRB( IVA) ,AVIRB( IVA+1) , 

1  (TEMP( I VC) , ISIG( I VC) , IVC=1 ,L) 

9060  FORMAT( 18X.2E12. 4,6(E12. 2,13)) 

6010  CONTINUE 

6020  CONTINUE 

6030  CONTINUE 
C 

END  IF 
C 

*1 .OUTPUT. 1383 

IF ( ICKV  .EQ.  1)  THEN 
WKILL  =  CZERO 
DO  112  J=1,LZMAX 
DO  111  IJ=1,100 

CHRNUM ( J , I J ) =C  ZERO 
CHRWAT( J , I J)=CZERO 
CHRENG(J,IJ)=CZERO 
CHRABS(J,IJ)=CZERO 
CHROUT( J, IJ)=CZERO 
CHRENT ( J , I J ) =C  ZERO 
111  CONTINUE 

112  CONTINUE 
ENDIF 

C 

IF( IVIRS  .EQ.  1)  THEN 
DO  115  IR=1 ,NRVIR 
DO  114  IC=1 ,NCVIR 
DO  113  IA=1 ,NAVIR 

SCOVOL( IA , IC , IR)=CZER0 

113  CONTINUE 

114  CONTINUE 

115  CONTINUE 
C 

ENDIF 

C  ■  v.  • 

*1, INPUT. 1853 

9051  WRITE(*,9070) 

9070  FORMAT 

1  (’OPREMATURE  EOF  WHILE  READING  CKV  OR  VIRS  DATA') 

GO  TO  9999 

9052  WRITE(*,9071) 

9071  FORMAT 

1  ( ' OILLEGAL  DATA  IN  FIELD  WHILE  READING  CKV  OR  VIRS  DATA') 

GO  TO  9999 
C 

*1 ,EHIST. 522 

IF  (ICKV  .EQ.  1)  THEN 

C  ***  SAVE  ZONE  INDEX  FOR  POSSIBLE  CALL  TO  TRANSRAD 
IRSAV  =  IR 
IRPSAV  =  IRPRIM 
DISADD=DIST0 

IFCMARK  .  NE.  1)  DISADD=DIST  +  l.D-07 
XCHR=X  +  DISADD*STH( 1)*CPH( 1) 

YCHR=Y  +  DISADD*STH(1)*SPH(1) 

ZCHR=Z  +  DISADD*CTH( 1 ) 

CALL  CKVRAD 
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oooooooooooooooooonooooono 


1  ( XCHR , YCHR , ZCHR , CTH , STH , CPH , SPH , W , T) 

C  ***  REINSTATE  ZONE  INDEX  FOR  POSSIBLE  CALL  TO  TRANSRAD 
IR  =  IRSAV 
IRPRIM  =  IRFSAV 
IF( IRPRIM  .LT.  0)  RETURN 


END  IF 
*1 ,DISTA. 133 

NBOSAV  =  MA(LSAVE) 

*D,PAREM.  6 

$  KLOOP,  LOOP,  I TYPE ,  NBOSAV 

* I, KNOCK. 367 

C  *CERENKOV  PRODUCTION  SELECTED  RANDOMLY  ALONG  A  SUBSTEP* 

0  A****************************************************** 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  CKVRAD 

1  ( XS , YS , ZS , CPOL , SPOL , CAZI , SAZI , WBETA , EBETA) 


VARIABLE  DESCRIPTIONS 


ACCEPT 

ALPHA 

ANGTCM 

C(l) 

C(  2) 

C(  3) 
CAFT 
CAZI 


FRACTION  OF  CERENKOV  RADIATION  ACCEPTED 
ZONE 

1. / 137. 0  (FINE  STRUCTURE  CONSTANT) 
l.E-8  (AN  ANGSTROM  IN  CENTIMETERS) 

X  PRESENT  DIRECTION  COSINE 
Y  PRESENT  DIRECTION  COSINE 
Z  PRESENT  DIRECTION  COSINE 

COSINE  OF  AZIMUTHAL  SCATTERING  ANGLE 


INTO 


C**********ceRENKOV  INPUTS******* 


CHER( 1 , J) 


MATERIAL  NUMBER  FOR  ZONE  J 

1  =  VACUUM 

2  =  AIR 


CHER( 2 , J) 
CHER( 3 , J) 

CHER(4 , J) 

CHER( 5 , J) 


*****CERENKOV 


3  =  C02 

4  =  GLASS 

5  =  QUARTZ 

EXTINCTION  COEFFICIENT  (1/CM) 

FRACTION  OF  LIGHT  TO  BE  MIRROR  REFLECTED 
FROM  ZONE  WALLS 

FRACTION  OF  LIGHT  TO  BE  DIFFUSE  REFLECTED 
FROM  ZONE  WALLS 

FRACTION  OF  LIGHT  TO  BE  ACCEPTED  FOR 
ENTRANCE  INTO  ZONE 
OUTPUTS****** 


CHRNUM(J,K) 


CHRWAT(J.K) 


CHRENG( J ,K) 


CHRABS( J , K) 


NUMBER  OF  CERENKOV  PHOTONS  PRODUCED  IN  ZONE  J 
FROM  WAVELENGTH  BIN  K 
PER  INCIDENT  SOURCE  PARTICLE 
CERENKOV  WEIGHT  PRODUCED  IN  ZONE  J 
FROM  WAVELENGTH  BIN  K  PER  INCIDENT 
SOURCE  PARTICE  WEIGHT  OF  1.  0 
CERENKOV  ENERGY  PRODUCED  IN  ZONE  J 
FROM  WAVELENGTH  BIN  K  PER  INCIDENT 
SOURCE  PARTICLE 
=  CHRENG(I)  +  EP*WP 

CERENKOV  WEIGHT  ABSORBED  (KILLED)  IN  ZONE  J 


oonoooooooo 


FROM  WAVELENGTH  BIN  K  PER 
INCIDENT  SOURCE  PARTICLE 
=  CHRABS(J)  +  WNEW  OR  WABS 
CHROUT(J.K)  CERENKOV  WEIGHT  THAT  EXITS  ZONE  J 

FROM  WAVELENGTH  BIN  K  PER 
INCIDENT  SOURCE  PARTICLE 
=  CHROUT(J)  +  WP 

CHRENT( J,K)  CERENKOV  WEIGHT  THAT  ENTERS  ZONE  (J+l) 

FROM  WAVELENGTH  BIN  K  PER 
INCIDENT  SOURCE  PARTICLE 
=  CHRENT( J+l)  +  WP 

C **************** *********** 


C  CN(J) 

C  CP(J) 

C  CP2 
C  POL 
C  CPRE 
C  CPRO 
C  CRAT 
C  CSQ 
C 

C  CX 
C  CZ 
C  DIFFU 
C  EBETA 
C  EE 
C  EP 
C  FST 
C 
C 

C  IBFIND 
C  J 

C  KALTOT 
C  KLIDE 
C  NR 
C  NRC 
C  NRPC 
C  NRMAX 
C  NZON 
C  PI 
C  PLANKC 
C  RA 

C  RANF(ARG) 

C 

C  RVIR 
C  SAFT 
C  SAZI 

C  SCOVOL( J, J, J) 
C 

C  SHD 
C  SIG 
C  SPECU 
C  SPOL 
C  ST 
C  STP 
C  V 


NORMAL  FROM  THIS  REGION  TOWARDS  NEXT  REGION 
DIRECTION  COSINES  IN  MOVING  FRAME 
COSINE  SQUARED 

COSINE  OF  POLAR  SCATTERING  ANGLE 
SUM  OF  SQUARES  OF  CURRENT  DIRECTION  COSINES 
PRODUCTION  COSINE  IN  THE  MOVING  FRAME 
RATIO  OF  COSINES 

SQUARE  ROOT  OF  SUM  OF  SQUARES  OF  X  AND  Y  DIRECTION 
COSINES 

X  DIRECTION  COSINE 
Z  DIRECTION  COSINE 
PROBABILITY  FOR  DIFFUSE  REFLECTION 
ELECTRON  ENERGY 
GAMMA  SQUARED 
ENERGY  OF  CERENKOV  PHOTON 
FRACTION  OF  LIGHT  TO  BE  TRANSMITTED  OUT 
THROUGH  ZONE  BOUNDARY 
■  1.0  -  CHER( 3 , J)  -  CHER(4 , J) 

FUNCTION  TO  FIND  NEXT  ZONE 
COUNTER 

MAXIMUM  NUMBER  OF  ALLOWED  REFLECTIONS 

COUNTER  FOR  NUMBER  OF  REFLECTIONS 

PRESENT  USER  REGION 

PRESENT  CODE  REGION 

NEXT  CODE  REGION 

MAX  ZONE  NUMBER 

ZONE  NUMBER 

3. 14159 

1. 2405E-10 

RANDOM  NUMBER 

RANDOM  NUMBER  PRODUCING  FUNCTION 

ARG  =  0  PRODUCES  UNIFORM  RANDOM  NUMBER 
VIRTUAL  SOURCE  RADIUS 

SINE  OF  AZIMUTHAL  SCATPERING  ANGLE 
DATA  FOR  VIRTUAL  SOURCE 
RADIUS  AND  TWO  ANGLES 
LENGTH  OF  ELECTRON  SUBSTEP  IN  CENTIMETERS 
PRODUCTION  OF  CERENKOV  PHOTONS  PER  CM  PER  MEV 
PROBABILITY  FOR  SPECULAR  REFLECTION 
SINE  OF  POLAR  SCATTERING  ANGLE 
DISTANCE  TO  BOUNDAPv 

SAMPLED  FRACTION  OF  ..LECTRON  SUBSTEP  IN  CM 
RATIO  OF  ELECTON  VELOCITY  TO  SPEED  OF  LIGHT  (BETA) 
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C  V2 
C  V2N2 
C  VARIND(I.J) 
C 

C  VOLVIR 
C  WABS 
C  WVLNTH(I) 

C 

C 

C 

C  WBETA 
C  WNEW 
C  VP 
C  WT(J) 

C  WZERO 
C  XB( J) 

C  XCHR(J) 

C  XPI(J) 

C  XPO(J) 

C  XREF 
C  XS 
C  XVIR 
C  YREF 
C  YS 
C  YVIR 
C  2E 
C  ZS 
C 


BETA  SQUARED 

BETA  SQUARED  TIMES  INDEX  OF  REFRACTION  SQUARED 
INDEX  OF  REFRACTION  FOR  MATERIAL  I  AND  WAVELENGTH 
BIN  J 

VOLUME  FOR  VIRTUAL  SOURCE 
CERENKOV  WEIGHT  ABSORBED  IN  ZONE 
WAVELENGTH  BIN  I 

THESE  BINS  RUN  FROM  1800  TO  5000  ANGSTROMS 
IN  200  ANGSTROM  WIDE  BINS  AND  THEN  FROM 
5500  TO  10000  ANGSTROMS  IN  500  ANGSTROM  BINS 
ELECTRON  WEIGHT 

REDUCED  WEIGHT  OF  CERENKOV  PHOTON 
WEIGHT  OF  CERENKOV  PHOTON 
TEMPORARY  POSITION  ARRAY 
INITIAL  WEIGHT  OF  CERENKOV  PHOTON 
TEMPORARY  POSITION  ARRAY 
PRESENT  POSITION 

POSITION  AT  BOUNDARY  JUST  INSIDE  THIS  REGION 
POSITION  AT  BOUNDARY  JUST  INSIDE  NEXT  REGION 
INDEX  OF  REFRACTION 

CURRENT  X  POSITION  (INPUT  IN  PARAMETER  LIST) 
VIRTUAL  SOURCE  X  POSITION 
INDEX  OF  REFRACTION  FOR  NEXT  ZONE 
CURRENT  Y  POSITION  (INPUT  IN  PARAMETER  LIST) 
VIRTUAL  SOURCE  Y  POSITION 
CHARGE  ON  SOURCE  PARTICLE  IN  ELECTRON  VOLTS 
CURRENT  Z  POSITION  (INPUT 
IN  PARAMETER  LIST) 


C********************************************** *********************** 


C 

*CA  CNSTNT 
*CA  PARAMS 
*CA  OUT 
*CA  CALC 
*CA  PAREM 
*CA  STTS 
*CA  GOMLOC 
*CA  DBG 
*CA  ORGI 
*CA  COG 
*CA  SECZN 

DOUBLE  PRECISION  DSEED ,  IRAN 
COMMON  /VAXRAN/  IRAN 

DIMENSION  C( 3) ,CN( 3) ,CP( 3) ,XCHR( 3) ,XPI ( 3) ,XPO( 3) 
DIMENSION  XCHR1(3) ,XB1(3) ,WT1(3) 

DIMENSION  WP(100) 

DATA  PLANKC , ANGTCM / 1 .  2405D-10.1.  D-8/ 

C  RANDOM  NUMBER  STATEMENT  FUNCTION  FOR  IBM 
RANF( DSEED )=  GGUBFS(DSEED) 

NCALL  =  NCALL  +  1 
ALPHA  =  1. DO/137. DO 
PI=DACOS( - 1. DO) 

TIME=1.  DO 
M=1 

C  CALCULATE  DIRECTION  COSINES  FOR  SOURCE  PARTICLE 
C( 1)  =  SPOL*CAZI 


O  D 


c 

c 


c 


c 


c 


100 


C( 2)  =  SPOL*SAZI 
C(3)  =  CPOL 

SAMPLE  RANDOM  LOCATION  FOR  CERENKOV  PRODUCTION  ALONG  A  SUBSTEP 
STP  =  SHD*RANF ( IRAN ) 

SET  ELECTON  LOCATION 
XCHR(l)  =  XS 
XCHR( 2)  =  YS 
XCHR(3)  =  ZS 

PRESENT  ZONE  OF  ELECTRON 
NRC  =  LBCZ 
NR  =  IROR(NRC) 

NRMAX  =  NZON 

UPDATE  RANDOMLY  SAMPLED  PRODUCTION  POSITION 
DO  100  J=1 , 3 
WT(J)  =  C(J) 

XB(J)  =  XCHR(J)  +  C( J)*(STP  -  SHD) 

XCHR(J)  =  XB(J) 

SAVE  VALUES  TO  INITIALIZE  MAIN  LOOP 
WT1(J)  =  WT(J) 

XBl(J)  =  XB( J) 

XCHRl(J)  =  XCHR(J) 

CONTINUE 

FIND  ZONE  OF  PRODUCTION  POSITION  AND  COMPUTE  UNCOLLIDED 
ESCAPE  DISTANCE. 

CALL  ZONEA 


C 

C  NEXT  ZONE  TO  BE  ENTERED.  IRPR.IM  AND  IR  SET  BY  SUBROUTINE  ZONEA 
NRC=IRPRIM 

C  PRESENT  ZONE  OF  PRODUCTION 


NR-IR 


C  SAVE  ZONE  TO  INITIALIZE  MAIN  LOOP 
NRC1  =  NRC 
NR1  =  NR 

C  INITIALIZE  XREF  TO  MINIMUM  VALUE 
MTRL=IDINT(CHER( 1 , NR ) ) 

XREF=VARIND( MTRL ,27) 

NTYPE=1 

C  ***  CALCULATE  VELOCITY  OF  SOURCE  PARTICLE  *** 

C  ELECTRON  ENERGY  NORMALIZED  TO  REST  MASS  ENERGY 
C  TOTAL  ENERGY  =  KINETIC  ENERGY  +  REST  MASS  ENERGY 
EE  =  1.  DO  +  EBETA/. 51 IDO 
C  NORMALIZED  VELOCITY  (V2  =  BETA  SQUARED) 

V2  =  1.D0  -  1.  D0/(EE*EE) 

C  (V  =  BETA) 

V  =  DSQRT(V2) 

C  (BETA  SQUARED  *  INDEX  OF  REFRACTION  SQUARED) 

V2N2  =  V2*XREF*XREF 

C  ***  DOES  PARTICLE  MEET  CERENKOV  THRESHOLD  *** 

IF  (XREF  .LE.  1.  DO)  GO  TO  1600 
IF(  V2N2  .LE.  1.  DO)  GO  TO  1600 
C  THE  CERENKOV  PHOTON  WEIGHT  TO  BE  USED  IS  THE 
C  WEIGHT  OVER  THE  RANGE  OF  WAVELENGTHS  FOR  THAT  BIN  AND  FOR  THE 
C  ENTIRE  LENGTH  OF  THE  SUBSTEP. 

DO  150  1=2 ,NMBINS 

WP( I-l)=2. D0*PI*ALPHA*SHD*WBETA*DABS((1. DO/ ( ANGTCM*WVLNTH( I - 1 ) ) 
1  -  1.  DO/(ANGTCM*WVLNTH( I)))* 
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n  o 


2  (1. DO-1.  DO/V/V/VARIND(MTRL, I -1)/VARIND(MTRL, 1-1) ) ) 

150  CONTINUE 
***  MAIN  LOOP 

SET  INITIAL  POSITION  AND  ZONE  OF  PRODUCTION 
DO  1560  K=1 ,NMBINS- 1 
DO  160  1=1,3 

XCHR(I)  =  XCHRl(I) 

WT( I)  =  WT1(I) 

XB( I )  =  XB1(I) 

160  CONTINUE 

NR  =  NR1 
NRC  =  NRC1 
IR  =  NR 
IRPRIM  =  NRC 

C  ***  GET  INDEX  OF  REFRACTION  *** 

200  XREF  =  VARIND(MTRL,K) 

C  ***  IS  INDEX  OF  REFRACTION  ABOVE  CUTOFF  FOR  CERENKOV?  *** 

IF(XREF  . LE.  1. DO)  GO  TO  1550 
C  CALCULATE  PHOTON  ENERGY  IN  MEV 

EP  =  PLAN'KC  /  ( ANGTCM*WVLNTH(  K  )  ) 

C  ***  CALCULATE  COSINE  FOR  CERENKOV  PHOTON  IN  MOVING  FRAME 
C  ***  =  1  /  (BETA  *  INDEX  OF  REFRACTION) 

C  ***  =  RATIO  OF  SPEED  OF  LIGHT  IN  MEDIUM  TO  ELECTRON  VELOCITY 

CPRO  =  l.D0/(V*XREF) 

C  ***  TALLY  PHOTON  NUMBER  COUNT,  WEIGHT,  AND  ENERGY  PRODUCED  HERE 
CHRNUM(NR,K)  =  CHRNUM(NR.K)  +  1.  DO 
CHRWAT(NR.K)  =  CHRWAT(NR,K)  +  WP(K) 

CHRENG(NR.K)  =  CHRENG(NR,K)  +  EP*WP(K) 

C  ***  BRING  POLAR  ANGLE  BACK  INTO  LAB  FRAME  *** 

CALL  FOLDC  CPOL , SPOL , CAZ I , SAZI , CPRO , CZ , SZ , CX , SX) 

C  ***  USE  SAME  INITIAL  PHOTON  WEIGHT 
WZERO  =  WP(K) 

C  ***  SET  DIRECTION  COSINES  AND  PREPARE  FOR  RAY  TRACE  *** 

C(  1 )  =  SZ-'CX 
C( 2)  =  SZ*SX 
C ( 3 )  =  CZ 

C  ***  CONDUCT  RAY  TRACE  OF  CERENKOV  PHOTON  THROUGH  GEOMETRY  *** 

NUMCT  =  NUMCT  +  1 
C  SET  INITIAL  COLLISITON  COUNTER 
KLIDE  =  0 

700  CALL  PATH( NR , NRC , XCHR , C , ST , NRPC , NSC , XPI , XPO , CN ) 

C  ***  CHECK  FORMATTING  HERE  *** 

C  ***  IF  NO  ZONE  IS  FOUND  THEN  PHOTON  ESCAPES 
IF( IRPRIM  . LT.  0)  THEN 
M  =  5 

GO  TO  1550 
ENDIF 
M  =  3 

C  ***  COMPUTE  THE  PHOTON  WEIGHT  ATTENUATION  AND  DETERMINE  THE  WEIGHT 
C  ***  LEFT  FOR  FURTHER  RAY  TRACE. 

WNEW  =  WP( K)*DEXP( -ST*CHER( 2 ,NR) ) 

C  ***  UPDATE  THE  WEIGHT  (CALCULATE  THE  WEIGHT  ABSORBED  IN  ZONE) 

WABS  =  WP(K)  -  WNEW 
VP(K)  =  WNEW 

C  ***  TALLY  THE  ABSORBED  WEIGHT  FOR  THAT  ZONE 
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CHRABSf NR,K)  =  CHRABS(NR.K)  +  WABS 
C  THIS  SECTION  CREATES  TALLIES  FOR  VIRTUAL  SOURCE 
IF( IVIRS  .EQ.  1)  THEN 

ICVIR  =  IBFIND(C( 3) ,NCVIR,CVIRB) 

IF( ICVIR  .  LE.  0)  GO  TO  1000 
SVIR  =  (HVIR  -  XCHR( 3))/C(3) 

XVIR  =  XCHR(l)  +  SVIR*C( 1) 

YVIR  =  XCHR( 2)  +  SVIR*C( 2) 

RVIR  =  DS  QRT (XVI R*XV I R  +  YVIR*YVIR) 

1RVIR  =  IBFIND(RVIR,NRVIR ,RVIRB) 

IF( IRVIR  .LE.  0)  GO  TO  1000 
CTVIR  =  1.  DO 
STVIR  =  0.  DO 
IF(RVIR  .GT.  0. DO)  THEN 
CTVIR  =  XVIR/RVIR 
STVIR  =  YVIR/RVIR 
ENDIF 

AVIR  =  0.  DO 

CSQ  =  DSQRT(C( 1)*C( 1)  +  C(2)*C(2)) 

IF(  CSQ  .LE.  0.  DO)  GO  TO  900 

AVIR  =  (CTVIR*C( 1)  +  STVIR*C( 2) )/CSQ 

IF(DABS(AVIR)  . GT.  1. DO)  AVIR  =  AVIR/DABS( AVIR) 

AVIR  =  DACOS(AVIR) 

900  IAVIR  =  IBFIND( AVIR,NAVIR,AVIRB) 

IF( IAVIR  .LE.  0)  GO  TO  1000 

VOLVIR  =  2.  DO*PI*(RVIRB( IRVIR  +  1)**2-RVIRB( IRVIR)**2) 

1  *(CVIRB( ICVIR  +  1)  -  CVIRBC ICVIR) ) 

2  *(AVIRB( IAVIR  +  1)  -  AVIRB( IAVIR) ) 

SCOVOL( IAVIR , ICVIR, IRVIR)= 

1  SCOVOL( IAVIR, ICVIR, IRVIR)+WP(K)/VOLVIR 

ENDIF 

C************  end  OF  VIRTUAL  SOURCE  SECTION  ***************** 

C  ***  INCREMENT  THE  NUMBER  OF  REFLECTIONS  (COLLISIONS  OF  PHOTON j 
1000  KLIDE  =  KLIDE  +  1 

C  ***  PUT  AN  UPPER  LIMIT  ON  THE  NUMBER  OF  REFLECTIONS  ALLOWED 
C  ***  (IE.  DON’T  GET  IN  AN  INFINITE  LOOP  BETWEEN  WALLS) 

IF( (KLIDE  .GT.  KALTOT).  OR. (WP(K).  LT. WZERO*l.  D-3) )  THEN 
WKILL  =  WKILL  +  WP(K) 

GO  TO  1550 
ENDIF 

C  ***  INITIALIZE  SPECULAR  AND  DIFFUSE  REFLECTION  COEFFICIENTS 
SPECU  =  0.D0 
DIFFU  =  0.  DO 

C  ***  SET  NEXT  REGION  INDEX 
NRP  =  NRPC 

C  ***  CHECK  TO  SEE  IF  WE  HAVE  A  NEXT  REGION 

C  ***  IF  SO  THEN  GET  THE  TRUE  SPECULAR  AND  DIFFUSE  COEFFICIENTS 
IF( NRPC  .GT.  0)  NRP=IROR( NRPC ) 

IF(NRP  .GT.  0)  THEN 
SPECU  =  CHER( 3 ,NRP) 

DIFFU  =  CHER(4,NRP) 

ENDIF 

C  ***  CHECK  TO  SEE  IF  PHOTON  UNDERGOES  SPECULAR  OR  DIFFUSE  REFLECTION 
RA  =  RANF( IRAN) -SPECU 
IF(RA  . LE.  0. DO)  THEN 
C  ***  SPECULAR  REFLECTION  *** 
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CPRE  =  C(1)*CN(1)  +  C( 2)*CN( 2)  +  C(3)*CN(3) 

1100  PS I  =  2. D0*CPRE 

C  ***  NEW  DIRECTION  =  OLD  DIRECTION  -  2*( NORMAL  COMPONENT) 

C  ***  VECTR2  =  VECTR1  -  2*(NMLVECT  DOT  VECTR1)*NMLVECT 
C(l)  =  C(l)  -  PSI*CN(1) 

C ( 2 )  =  C( 2)  -  PSI*CN( 2) 

C( 3)  =  C( 3)  -  PSI*CN(3) 

C  ***  SET  NEW  POSITION 

DO  1200  J=1 , 3 

XCHR(J)  =  XPI(J) 

1200  CONTINUE 

C  ***  REPEAT  UNTIL  PHOTON  TERMINATES  *** 

GO  TO  700 
END  IF 

C  ***  CHECK  FOR  DIFFUSE  REFLECTION  *** 

RA  =  RA  -  DIFFU 
IF(RA  . LE.  0. DO)  THEN 

C  ***  SAMPLE  RANDOM  POLAR  COSINE  RELATIVE  TO  SURFACE  NORMAL 
CP( 3 )  =  - RANF ( I RAN ) 

C  ***  SINE  OF  POLAR  SCATTERING  ANGLE 

SINE  =  DSQRT( 1. DO  -  CP(3)*CP(3)) 

C  ***AZIMUTHAL  SCATTERING  ANGLE  RANDOM  BETWEEN  +  OR  -  PI 
THE  =  P I * ( 2 . D0*RANF ( I RAN )  -  1.  DO) 

CP(1)  =  S I NE*DCOS ( THE ) 

CP( 2)  =  S I NE*D S I N ( THE ) 

C  ***  ROTATE  DIRECTION  COSINES  BACK  INTO  LAB  FRAME 
CALL  ROTCPC( CN , CP , C ) 

C  ***  UPDATE  CURRENT  POSITON 
DO  1300  J=1 , 3 

XCHR(J)  =  XPI(J) 

1300  CONTINUE 

GO  TO  700 
END  IF 

C  ***  TRANSMIT  THE  CERENKOV  PHOTON  *** 

C  ***  SET  CUSv'  OF  INCOMING  RAY  RELATIVE  TO  SURFACE  NORMAL 
CPRE  =  C( 1)*CN( 1)  +  C(2)*CN( 2)  +C(3)*CN(3) 

YREF  =  1.  DO 
C  ***  SET  ZONE  NUMBER 
NRP  =  NR PC 

IFCNRPC  . GT.  0)  NRP=IROR( NRPC ) 

ACCEPT  =  1. DO 

C  ***  IF  THERE  IS  A  NEXT  ZONE  SET  INDEX  OF  REFRACTION  AND  FRACTION 
C  ***  ACCEPTED 

IF( NRP  .GT.  0)  THEN 

MTRL1  =  IDINT( CHER( 1 ,NRP) ) 

YREF  =  VARIND( MTRL1 ,K) 

ACCEPT  =  CHER( 5 ,NRP) 

END  IF 

C  ***  TEST  TO  SEE  IF  THE  PHOTON  IS  REFLECTED  OR  TRANSMITTED 
IF(RANF( IRAN)  . GT.  ACCEPT)  GO  TO  1100 
C  ***  DOES  SNELLS  LAW  ALLOW  TRANSITION?  (SNELL’S  LAW  IN  TERMS  OF 
C  ***  COSINE  OF  ANGLE  SQUARED) 

CP2  =  1.  DO  -  XREF*XREF*( 1 . DO  -  CPRE*CPRE)/( YREF* YREF) 

C  ***  IF  COSINE  OF  ANGLE  SQUARED  LESS  THAN  OR  EQUAL  ZERO  THEN 
C  ***  SPECULAR  REFLECTION  OCCURS  (OUTSIDE  CRITICAL  ANGLE; 

IF( CP2  . LE.  0. DO)  GO  TO  1100 
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SPRE  =  DSQRT( 1. DO  -  CPRE*CPRE) 

CAFT  =  DSORT(CP2) 

SAFT  =  DSCjRT(  1.  DO  -  CP2) 

SRAT  =  SAFT/ SPRE 
CRAT  =  CAFT  -  SRAT*CPRE 

C  ***  DIRECTION  COSINES  AFTER  ENTERING  NEW  ZONE  IN  LAB  FRAME 
DO  1400  J=1 . 3 

C(J)  =  SRAT*C(J)  +  CRAT*CN(J) 

1400  CONTINUE 

C  ***  INCREMENT  CERENKOV  WEIGHT  EXITING  PREVIOUS  ZONE 
CHROUT(NR,K)  =  CHROUT(NR.K)  +  WP(K) 

C  ***  is  THERE  A  NEXT  ZONE? 

IF(NRP  .GE.  0)  THEN 

C  ***INCR£M£NT  CERENKOV  WEIGHT  ENTERING  NEXT  ZONE 
CHRENT(NRP ,K)  =  CHRENT(NRP ,K)  +  WP(K) 

C  ***  UPDATE  ZONE  AND  INDEX  OF  REFRACTION 
NR  =  NRP 
NRC  =  NRPC 
XREF  =  YREF 
C  ***  UPDATE  POSITION 

DO  1500  J=1 , 3 

XCHR(J)  =  XPO(J) 

1500  CONTINUE 

C  ***  REPEAT  UNTIL  ALL  ZONES  COMPLETED 

IF(NRP  .  LT.  NRMAX)  GO  TO  700 
ENDIF 

1550  CONTINUE 
1560  CONTINUE 
1600  RETURN 
END 

C  ***  SUBROUTINE  TO  ROTATE  COORDINATES  BACK  INTO  LAB  FRAME 
SUBROUTINE  ROTCPC(CP,D,C) 

IMPLICIT  REAL* 8  (A-H.O-Z) 

DIMENSION  CP( 3) ,D( 3) ,C( 3) ,E3( 3) ,E2( 3) ,E1( 3) 

DO  100  J=1 , 3 
E3( J)  =  CP(J) 

100  CONTINUE 

E22  =  E3( 1)*E3( 1)  +  E3(2)*E3(2) 

IF( E22  .  LE.  0. DO)  THEN 
E 2  ( 1 )  =  1.  DO 
E  2 ( 2)  =  0. DO 
ELSE 

E22  =  DSQRT(E22) 

E2( 1)  =  -E3(2)/E22 
E2( 2)  =  E3(l)/E22 
ENDIF 

E  2 ( 3 )  =  0. DO 

E 1 ( 1 )  =  E 2 ( 2)*E3( 3)  -  E2(3)*E3(2) 

El( 2)  =  E2( 3)*E3( 1)  -  E2(1)*E3(3) 

El( 3)  =  E2(1)*E3(2)  -  E2(2)*E3(1) 

DO  200  J=1 , 3 

C(J)  =  D(1)*E1(J)  +  DC  2)*E2( J)  +  D( 3 )*E3( J) 

200  CONTINUE 
RETURN- 
END 
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SUBROUTINE  PATH( NR , NRG , X , C , ST , NRGP , NSCP , XPI , XPO , CN) 

C  ***  RAY  TRACING  SUBROUTINE  **** 

*CA  CNSTNT 
*CA  PARAMS 
*CA  PAREM 
*CA  STTS 
*CA  GOMLOC 
*CA  DBG 
*CA  ORGI 
*CA  SEC2N 
*CA  COG 

COMMON/NRML/UN(  3) 

DIMENSION  X( 3) ,C( 3) ,XPO( 3) ,XPI( 3) ,CN( 3) 

C  ***  ACCEPT  VERSION 
DIST  =  CZERO 
DISTO  =  PINT 
IRPRIM  =  NRG 
IR  =  NR 

KLOOP  =  KLOOP  +  1 
DO  100  1=1,3 
XB( I )=X( I ) 

WT( I)=C(I) 

100  CONTINUE 
ICALL  =  0 
IZE  =  0 
CALL  DISTA(S) 

IF( IRPRIM  .GE.  0)  THEN 
ST  =  DIST 
DO  200  1=1,3 

XPI(I)  =  X(I)  +  (ST*0. 99999D0)*C( I) 

XPO(I)  =  X(:)  +  (ST*1.00001D0)*C(I) 

200  CONTINUE 

NRGP  =  IRPRIM 

NSCP  =  LSURF 

CALL  NORML( NRG , LSAVE ) 

C  ***  SET  NORMAL  TO  DIRECTION  UNTIL  NORMAL  IS  AVAILABLE 
DO  300  1=1,3 

CN( I )  =  -UN(I) 

300  CONTINUE 

END  IF 
RETURN 
END 

SUBROUTINE  NORML( IRR ,NASC) 

C  ***  USER  BEWARE:  MUCH  OF  THIS  ROUTINE  IS  UNTESTED.  IT  IS  ASSUMED 
C  ***  THAT  XB+WT*DIST  IS  ON  SURFACE  LSURF  OF  BODY  NASC.  THUS  A  NORMAL 
C  ***  CALL  MUST  BE  PROCEEDED  BY  A  CALL  TO  G1  OR  LOOKZ. 

*CA  CNSTNT 
*CA  PARAMS 
*CA  PAREM 
*CA  GOMLOC 
*CA  DBG 
*CA  ORGI 
*CA  SECZN 
*CA  COG 

COMMON  /NRML/UN(  3) 

DIMENSION  XP( 3) ,WP(3) ,X(3) ,H( 3) , IB0X1( 3) , IB0X2( 3) 
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DATA  IBOXl/10,7,4/ 

DATA  IB0X2/4, 10,7/ 

LSUR  =  IABS(LSURF) 

DO  100  1=1,3 

XP(I)  =  XB( I )  +  DIST*WT(I) 

100  CONTINUE 

N  =  LOCREG(IRR) 

NUM  =  NUMB0D(IRR)*4  +  N  -  4 
DO  200  I=N ,NUM ,4 

IF(NASC  .EQ.  I)  GO  TO  400 
200  CONTINUE 

WRITE( 6 , 300)  IRR ,NASC ,  (MA( I ) , I=N,NUM,4) 

300  FORMATC IX,' INVALID  REGION  OR  BODY  IN  NORMAL.  IR  ='JI7,,NASC 
1  17,/, (1017)) 

CALL  PR( 1) 

STOP 

400  L  =  MA( 1+1 ) 

I TYPE  =  MA(L+1) 

K  =  MA( L+4) 

C  ARB  SPH  RCC  REC  TRC  ELL  BOX  WED  RPP 

GO  TO  (500,1200,1400,1500,1400,1600,1700,1500,1900),  I TYPE 
C  ***  ARB 

500  I=K  +  1  +  4*( LSUR  -  1) 

DO  1110  J=1 , 3 

UN(J)  =  FPD( I+J) 

1110  CONTINUE 

GO  TO  2000 

C  ***  SPH 

1200  DO  1210  1=1,3 

UN( I )  =  (XP(I)  -  FPD( K+I+l) )/FPD( K+5) 

1210  CONTINUE 

GO  TO  2000 
C  ***  TRC  AND  RCC 
1400  H2=0. DO 

DO  1410  1=1,3 

H( I )=FPD( K+4+I ) 

H2=H2  +  H(I)*H(I) 

1410  CONTINUE 

GO  TO  (1415,1415,1435),  LSUR 
1415  H2  =  DSQRT(H2) 

DO  1420  1=1,3 

UN( I )  =  H(I)/H2 
1420  CONTINUE 

GO  TO  2000 
1435  DO  1440  1=1,3 

X(I)  =  XP(I)  -  FPD(K+1+I) 

1440  CONTINUE 

XMU  =  0.  DO 
A2  =  0.  DO 
DO  1445  1=1,3 

J  =  MOD( 1,3)  +  1 

A2  =  A2  +  (H(I)*X(J)  -  H(J)*X(I))**2 
XMU  =  XMU  +  FPD( K+4+I )*X( I ) 

1445  CONTINUE 

R  =  FPDfK+8)  -  FPD( K+9 ) 

IF( ITYPE  .EQ.  3)  R  =  0.  DO 
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1450 

1500 

1501 

C  *** 
1600 

1620 

1630 

1635 

1650 

1651 

C  *** 

1700 

1720 

1730 

1735 


B1  =  R/DSQRT(H2*(H2  +  R*R) ) 

B2  =  1. D0/DSQRT(A2*(H2  +  R*R) ) 

DO  1450  1=1,3 

UN( I)  =  B2*(X( I )*H2  -  H(I)*XMU)  +  B1*H(I) 

CONTINUE 

GO  TO  2000 

WRITE( 6 , 1501)  I TYPE 

F0RMAT(1X, 'NOT  IMPLEMENTED  IN  NORMAL.  I TYPE  =',I5) 

STOP 

ELL 

DX  =  0.  DO 
AN  =  0.  DO 

C2  =  FPD(K+8)**2/4.  DO 
DO  1620  1=1,3 
IA  =  I  +  K  +  1 

X(I)  =  XP(I)  -  . 5D0*(FPD(IA)  +  FPD( IA+3) ) 

H( I )  =  .  5D0*(FPD(IA+3)  -  FPD(IA)) 

DX  =  DX  +  X( I )*H( I ) 

CONTINUE 
DO  1630  1=1,3 

UN( I)  =  C2*X( I)  -  DX*H( I ) 

AN  =  AN  +  UN( I )*UN( I ) 

CONTINUE 

IF( AN  .GE.  l.D-12)  THEN 
AN  =  DSQRT(AN) 

DO  1635  1=1,3 

UN( I)  =  UN( I) /AN 
CONTINUE 
GO  TO  2000 
END  IF 

WRITEC 6 , 165 1 )  NBO,LSURF,(XP(I) ,1=1,3) 

FORMAT( IX, 'ROUND-OFF  ERROR  IN  NORMAL.  NBO=' , 15 , ' LSURF=’ , 15 , 
L  'XP=' ,3E12.  5) 

STOP 

BOX 

L  =  1  +  ( LSUR  -  l)/2 

11  =  IBOXl(L) 

12  =  IB0X2(L) 

DO  1720  1=1,3 

X(I)  =  FPD(K+I1+I) 

H( I )  =  FPD( K+I2+I ) 

CONTINUE 
AN  =  0. DO 
DO  1730  1=1,3 

J  =  MOD( 1,3)  +  1 

JJ  =  MOD( 1+1,3)  +  1 

UN( I )  =  X( J)*H( JJ)  -  X( JJ)*H( J) 

AN  =  AN  +  UN(I)*UN(I) 

CONTINUE 

IF( AN  .LT.  1.0D-12)  GO  TO  1650 
AN  =  DSQRT(AN) 

DO  1735  1=1,3 

UN ( I )  =  UN(I) /AN 
CONTINUE 
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GO  TO  2000 
C  ***  RPP 

1900  DO  1905  1=1,3 

UN( I)  =  0.  DO 
1905  CONTINUE 

GO  TO  (1910,1910,1920,1920,1930,1930),  LSUR 
1910  UN( 2)  =  1.  DO 
GO  TO  2000 
1920  UN( 1)  =  1. DO 
GO  TO  2000 
1930  UN( 3)  =  1.  DO 
2000  XMU  =  0.  DO 

DO  2005  1=1,3 

XMU  =  XMU  +  WT( I )*UN( I ) 

2005  CONTINUE 

XMU  =  - 1. DO*XMU 
C  AN  =  SIGN( 1. DO,  XMU) 

IF(XMU  .LT.  0.  DO)  THEN 
AN  =  - 1.  DO 
ELSE 

AN  =  1.  DO 
END  IF 

XMU  =  -l.DO*XMU 
DO  2010  1=1,3 

UN(I)  =  AN*UN( I ) 

2010  CONTINUE 
RETURN 
END 

FUNCTION  IBFIND(X ,MX,XB) 

C  ***  IBFIND  FINDS  CORRECT  TALLY  BINS  *** 
IMPLICIT  REAL* 8  (A-H.O-Z) 

DIMENSION  XB(IOO) 

IBFIND=0 

IF( X  . LT.  XB(1))  GO  TO  7777 
DO  1000  1=1, MX 

IF(X  .LE.  XB( 1+1) )  GOTO  1010 
1000  CONTINUE 

GO  TO  7777 
1010  IBFIND=I 
7777  RETURN 
END 
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APPENDIX  C.  SAMPLING  SUBROUTINE  LISTING 


* I DENT  DEFINE 
♦DEFINE  IBM 
♦DEFINE  DOUBLE 
♦DEFINE  ACCEPT 
♦I DENT, PHILLIPS 
♦IDENT.PERMT 
*D  P ARAMS.  17 

PARAMETER  (KPTMAX=10000 ,  INSTAT=100, 

♦I, HIST. 94 

COM  =  CZERO 
♦I, HIST. 223 

IF( ( IB  . EQ.  1)  .AND.  (IHIST  .  LE.  50))  THEN 
I F ( IHIST  .EQ.  1)  WRITE( 6,801) 

C  SEND  ADITIONAL  OUTPUT  TO  OUTPUT  UNIT  6 

801  FORMAT ( IX , 1HI , 2X , 2HLB , 7X , 1HX , 14X , 1HY , 14X , 1HZ , 

1  13X , 1HU , 13X , 1HV , 12X , 1HW, 12X, 3HCOM, 10X, 7HENERGY  /) 

WRITE(6 , 802)  IHIST, LB ,XB,WT, COM, T 
802  FORMAT (2I3,8(2X,1PE12. 5)) 

C 

END  IF 
C 

♦IDENT.CKVTR 
♦I, CALC.  78 

6  , AIRP,C02P,TMPKEL 

♦I, OUT.  30 

f  CERENKOV  ADDITIONS  TO  COMMON  BLOCK 
COMMON  /OUT/ 

1  CHER( 5 , 100) ,CHRNUM( 100 , 100) , 

2  CHRWAT( 100 , 100) ,CHRENG( 100, 100) ,CHRABS( 100 , 100) , 

3  CHROUTC 100 , 100) ,CHRENT( 100,100) ,TRNUM( 100,100) , 

4  TRWAT ( 100,100) ,TRENG( 100,100) ,TRABS( 100,100) , 

5  TROUT ( 100,100) , TRENT ( 100,100) ,NRVIR,NCVIR,NAVIR, 

6  CVIRB( 11) , AVIRB( 11) ,RVIRB( 11) ,VARIND( 100 , 100) ,WVLNTH( 100) , 

7  SCOVOL( 10,10,10) , HVIR , ITRAD ,NRBOTH , ICKV, IVIRS ,KALTOT, 

8  LBTRAD , NBOTRD , JSRTRD , NCALL , NMB INS , WKILL 
♦I, INPUT.  54 

DIMENSION  WVLNT( 100) ,VARIN( 100,100) 

iV***************************************************************** 


,,*****THiS  SECTION  ASSIGNS  VALUES  TO  THE  TABLE  OF  INDEX  OF  ***** 
* ♦♦♦♦♦REFRACTION  AND  WAVELENGTH.  THE  TABLES  CREATED  ARE  ***** 

C *****STRUCTURED  TO  COVER  WAVELENGTHS  IN  200  ANGSTROM  BINS  ***** 

*****FROM  1800  THROUGH  5000  ANGSTROMS,  AND  IN  500  ANGSTROM  ***** 

C  '*••-** BINS  FROM  5500  THROUGH  10000  ANGSTROMS.  ***** 

*****THE  VARIND  TABLr.  CONTAINS  INDICES  OF  REFRACTION  FOR  ***** 
♦♦♦♦♦EACH  MATERIAL.  ***** 

♦  'VcVrVrA’  VrVnWrVr 

***** V AR I ND ( I , J )  --  INDEX  OF  REFRACTION  FOR  MATERIAL***** 

*****  MATERIAL  I  ASSOCIATED  WITH  ***** 

*****  WAVELENGTH  BIN  J  ***** 

*  Vf  Vr  *V  Vr  Vr  Vr  Vr  Vr 

*****WVLNTH(J)  --  WAVELENGTH  VALUE  FOR  BIN  J  ***** 
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UAVFT  FNftTU^ 

DATA  (WVLNT( I), 1=1, 27)/  1800,2000, 

1  2200,2400,2600,2800,3000, 

2  3200,3400,3600,3800,4000, 

3  4200,4400,4600,4800,5000, 

4  5500,6000,6500,7000,7500, 

5  8000,8500,9000,9500,10000  / 

C  VACUUM 

DATA  (VARIN( 1,1) ,I=1,27)/27*1. DO/ 

C  AIR 

DATA  (VARIN(2, I) ,1=1,27)/ 

1  1. 00034D0 , 1. 00031D0, 1. 00030D0,1. 00030D0.1. 00029D0, 

2  1.  00029D0,1.  00029D0 , 1.  00029D0.1. 00028D0,1. 00028D0, 

3  1. 00028D0 , 1. 00028DO , 1. 00028D0,1. 0002800,1. 00028D0, 

4  1.  00028D0 , 1.  00028D0 , 1. 00028DO,1. 00028D0.1. 00028D0, 

5  1.  00028D0 , 1.  00027D0 , 1.  00027D0.1. 00027D0.1. 00027DO, 

6  1. 00027D0 , 1. 00027D0  / 

C  CARBON  DIOXIDE  (STANDARD  ATMOSPHERE) 

DATA  ( VARIN( 3, I) ,1=1,27)/ 

1  1.  00055D0 , 1.  0005200,1.  0005100,1. 00050D0.1. 00049D0, 

2  1. 00048D0 , 1. 00048D0 , 1. 0004700,1. 0004700,1. 00046D0, 

3  1.  00046D0 , 1.  0004600,1.  00046D0.1. 00046D0,1.  00045D0, 

4  1. 00045D0 , 1.  00045D0 , 1. 00045D0.1. 00045D0.1. 00045D0, 

5  1.  00045D0 , 1.  00045D0 , 1.  00045D0.1. 00044D0.1. 00044D0, 

6  1. 00044D0.1. 00044DO  / 

C  GLASS  TYPE  458/678(MIL-G-174)  (FUSED  SILICA  DENSITY  2.202  G/CM**3) 
DATA  (VARIN(4, I) ,1=1,27)/ 

1  1.  58529D0 , 1. 5505 IDO , 1. 52845D0,1. 51333D0.1. 50352D0, 

2  1. 49565D0, 1. 4877900,1. 48274D0,1. 47884D0.1.  47554D0, 

3  1. 47283D0 , 1. 47012D0.1. 46830D0,1. 46648D0, 1. 46492D0, 

4  1. 483f^00, 1. 46233D0 , 1. 45991D0, 1. 45804D0, 1. 45664D0, 

5  1.  4552  vL)0 , 1.  4542400,1.  45332D0.1.  45250D0,1.  45175D0, 

6  1. 45107D0.1. 45042D0  / 

C  CRYSTAL  QUARTZ 

DATA  (VARIN(5, I) ,1=1,27)/ 

1  1. 6600600,1.  65087D0,1. 62626D0, 1. 61011D0, 1. 6015800, 

2  1. 5930600,1. 5845300,1. 57600DO,1. 56747D0.1. 56413D0, 

3  1.  56080D0,1.  55779D0 , 1. 55554D0, 1. 55349D0, 1. 55194D0, 

4  1. 55039D0 , 1. 54884D0,1. 54616D0,1. 54393D0.1. 54247D0, 

5  1. 5410100,1.  53955D0 , 1.  53839D0.1. 53744D0,1.  53663D0, 

6  1. 5358100,1. 53502D0  / 

C 

*1, INPUT.  436 

C  CERENKOV  ADDITIONS  TO  DEFAULT  INPUT  PARAMETERS 
NMBINS=27 
NCALL=0 
ICKV=0 
IVIRS=0 
ITRAD=0 

C  END  OF  CERENKOV  ADDITIONS 
*1, INPUT. 575 

Q*****-V****-.V-V**tV******* 
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C  CERENKOV 

Q  irir  Vc  tV-jV  iric  Jr ?V*Ar*  ir'k •iHHck  *  *  Jr Jrfr 

ELSE  IF  (K0P( 'CERENKOV' )  .  GE.  0)  THEN 
I MATCH  =  1 
ICKV  =  1 

C  SET  PRESSURE  FOR  AIR  IN  STANDARD  ATMOSPHERES 
AIRP  =  PARM(l) 

C  SET  PRESSURE  FOR  C02  IN  STANDARD  ATMOSPHERES 
C02P  =  PARM( 2) 

C  SET  TEMPERATURE  IN  KELVIN,  THEN  SCALE  TO  293  KELVIN 
TMPKEL  =  PARM( 3) 

TMPKEL  =  TMPKEL/ 29 3.  DO 
DO  93774  1 1=1 ,NMBINS 

VARINDC 1 , 1 1 )=VARIN( 1,11) 

C  INDICES  OF  REFRACTION  -  1  CHANGE  PROPORTIONAL  TO  PRESSURE  CHANGE 
VARINDC 2 , Il)=( VARIN( 2 , 1 1 ) -1.  D0)*AIRP/TMPKEL  +  1. DO 
VARINDC  3 , 1 1 )=( VARIN( 3 , 1 1) - 1.  D0)*C02P/TMPKEL  +  1.  DO 
93774  CONTINUE 

DO  93884  Jl=4 , 20 

DO  93874  1 1=1 .NMBINS 

VARINDC J1,I1)=VARIN(J1, II) 

93874  CONTINUE 

93884  CONTINUE 

DO  93894  1 1=1 .NMBINS 

WVLNTH C 1 1 )  =WVLNT(  1 1 ) 

93894  CONTINUE 

KALTOT  ~  PARMC4) 

DO  3034  JJJ=1,N20N 

CALL  OPREADC 1, INUNIT, IECHO, IVAL.EOFLAG) 

C 

C  CHECK  INPUT  FILE  FOR  COMPLETE  OPTICAL  DATA 
C 

IFCEOFLAG)  THEN 

WRITEC*,*)  'PREMATURE  END  OF  OPTICAL  DATA' 

STOP 

EN'DIF 

IF  CKOPC '+' )  -GE.  0)  THEN 
C 

C  ASSIGN  VALUES  TO  CERENKOV  INPUT  MATERIAL  NUMBER  FOR  ZONE  JJJ 

CHERC 1 , JJJ)  =  PARMC1) 

C  EXTINCTION  COEFFICIENT  Cl/CM)  --  ATTENUATION  COEFFICIENT 

CHER( 2 , JJJ)  =  PARMC2) 

C  FRACTION  OF  LIGHT  TO  BE  MIRROR  REFLECTED  FROM  ZONE  WALLS 

CHERC 3, JJJ)  =  PARMC  3) 

C  FRACTION  OF  LIGHT  TO  BE  DIFFUSE  REFLECTED  FROM  ZONE  WALLS 

CHERC  4 , JJJ)  =  PARMC 4) 

C  FRACTION  OF  LIGHT  TO  BE  ACCEPTED  FOR  ENTRANCE  INTO  ZONE 

CHERC 5, JJJ)  =  PARMC 5) 

C 

ELSE 

WRITEC*, 8903) 

8903  FORMATC IX, ' ILLEGAL  CKV  OR  VIRS  ZONE  DATA') 

GO  TO  9999 
END  IF 

3034  CONTINUE 

Q  V?  *  iri f  V?  *********  -.V  *  *  *  *  *  * 


C  VIRTUAL- SOURCE 

C********************* 

ELSE  IF  (KOP( ’VIRTUAL-SOURCE' )  .  GE.  0)  THEN 
IMATCH  =  1 
IVIRS  =  1 
HVIR  =  PARM(l) 

NRVIR  =  PARMC  2) 

NCVIR  =  PARM( 3) 

NAVIR  =  PARM( 4) 

RE AD( INUNIT, *,END=9051,ERR=9052) 

1  (RVIRB(IR) ,IR=1, NRVIR)  ,RVIRB(NRVIR  +  1) 

READ( INUNIT, * ,END=9051 ,ERR=9052) 

1  (CVIRB(IC) ,IC=1, NCVIR), CVIRB(NCVIR  +  1) 

READ( INUNIT, *,END=9051 ,ERR=9052)(AVIRB( IA) , 

1  IA=1, NAVIR)  ,AVIRB( NAVIR  +  1) 

C ********************** 

C  TRANSITION  RADIATION 

Q***Vf******iV')V***,!Wf**iV*')V 

ELSE  IF(KOP( 'TRANSITION-RADIATION')  .  GE.  0)  THEN 
IMATCH  =  1 
I TRAD  =  1 
LBTRAD  =  PARMC 1) 

NBOTRD  =  PARM( 2) 

JSRTRD  =  PARMC 3) 

*1, INPUT. 1178 

IFCICKV  .EQ.  1)  THEN 

C  CERENKOV  OPTICAL  AND  VIRTUAL  SOURCE  DATA . 

WRITEC  6 ,9005 ) 

9005  FORMAT(6H  ZONE,6X,6H  MTRL  ,6X,6HABSORB,6X,6HMIRROR,5X,7HDIFFUSE 
1  , 6X , 6HACCEPT , 4X , 8HTRANSMIT) 

DO  90051  1=1 , LZMAX 

C  FRACTION  OF  LIGHT  TO  BE  TRANSMITTED  OUT  THROUGH  ZONE  BOUNDARY 
FST  =  1  -  CHER( 3,1)  -  CHERC4.I) 

WRITE ( 6, 90050)1, ( CHER( J , I ) , J=1 , 5) ,FST 

90050  FORMATC 16 , 1P6E12.  4) 

90051  CONTINUE 
C 

END  IF 
C 

IF( ITRAD  .EQ.  1)  THEN 

WRITEC  6,90311)  LBTRAD , NBOTRD , JSRTRD 
90311  FORMATC/'  TRANSITION  RADIATION  LOCATION,  ZONE  BODY  SURFACE' 

1  ,316) 

C 

ENDIF 

C 

IF  C IVIRS  .EQ.  1)  THEN 

WRITEC 6, 90300)  HVIR, NRVIR , NCVIR , NAVIR 
90300  FORMATC /17H  VIRTUAL  SOURCE  Z,  E12.4 

1  I'  VIRTUAL  SOURCE  TALLY  BINS,  RADIAL,  POLAR,  AZIMUTH,’ 

2  ,  316) 

WRITEC  6 ,9050) (RVIRBC IR) , IR=1 , NRVIR) ,RVIRB( NRVIR  +  1) 

WRITEC  6 , 9055 )CCVIRB( IC) , IC=1 , NCVIR) ,CVIRB( NCVIR  +  1) 

WRITEC  6 , 9060) ( AVIRBC IA) , IA=1 , NAVIR) ,AVIRB( NAVIR  +  1) 

9050  FORMATC 8H  RADIAL  .10E12.4) 

9055  FORMATC 8H  POLAR  , 10E12.4) 
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9060  FORMAT(8H  AZIMUTH  .10E12.4) 

C 

ENDIF 

C 

*I,PREP. 659 

IF( ICKV  .EQ.  1)  THEN 
WKILL  =  CZERO 
DO  241  LLL=1 , 100 
DO  2401  JJJ=1,100 

CHRNUM(LLL.JJJ)  =  CZERO 
CHRWAT(LLL.JJJ)  =  CZERO 
CHRENG(LLL.JJJ)  =  CZERO 
CHRABS(LLL.JJJ)  =  CZERO 
CHROUT(LLL, JJJ)  =  CZERO 
CHRENT(LLL.JJJ)  =  CZERO 
2401  CONTINUE 

241  CONTINUE 
C 

ENDIF 

C 

IF( IVIRS  .EQ.  1)  THEN 
DO  244  IR=1 ,NRVIR 
DO  243  IC=1 ,NCVIR 
DO  242  IA=1 ,NAVIR 

SCOVOL( IA, IC , IR)  =  CZERO 

242  CONTINUE 

243  CONTINUE 

244  CONTINUE 
C 

ENDIF 

C 

*1, OUTPUT.  47 

DIMENSION 

1  CHRN( 100) , CHRW ( 100) ,CHRE( 100) ,CHRA( 100) ,CHRO( 100) ,CHRI( 100) 
INTEGER  ISIGA( 100) ,ISIGB( 100) , ISIGCC 100) , ISIGD( 100) ,ISIGE( 100) , 

1  ISIGF( 100) , IWVL( 100) , ISIG1( 100) , ISIG2( 100) , ISIG3( 100) , 

2  ISIG4(100),ISIG5(100),ISIG6(100) 

*1, OUTPUT.  561 

IF( ICKV  .EQ.  1)  THEN 
WTESC  =  CZERO 
WTABS  =  CZERO 
VTOT  =  CZERO 

C  SCALE  RESULTS  TO  RESULTS  PER  HISTORY 
WKILL  =  WKILL/CIMAX 
DO  66  1=1 ,LZMAX 

DO  56110  J=1 .NMBINS- 1 

CHRNUM( I , J ) =( CHRNUM( I , J ) / ( C I MAX ) ) 

CHRWAT( I , J)=CHRWAT( I , J)/CIMAX 
CHRENG( I , J ) =CHRENG( I , J ) /C IMAX 
CHRABS( I , J)=CHRABS( I , J)/CIMAX 
CHROUT( I , J)=CHROUT( I , J)/CIMAX 
CHRENT( I , J)=CHRENT( I , J)/CIMAX 
56110  CONTINUE 

66  CONTINUE 

DO  10225  1=1,100 
CHRN( I )=CZERO 
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CHRW( I )=CZERO 
CHRE( I )=CZERO 
CHRA( I )=CZERO 
CHRO( I )=CZERO 
CHRI( I)=CZERO 
10225  CONTINUE 

C  COMPUTE  THE  TOTALS  FOR  EACH  ZONE 
DO  11230  1=1 , LZMAX 

DO  11225  J=1 ,NMBINS-1 

CHRN ( I ) =CHRN ( I ) +C HRNUM ( I , J ) 

CHRW( I )=CHRW( I )+CHRWAT( I , J) 

CHRE ( I ) =CHRE ( I ) +CHRENG( I , J ) 

CHRA( I)=CHRA( I)+CHRABS( I , J) 

CHRO( I)=CHRO( I)+CHROUT( I, J) 

CHRI( I )=CHRI( I)+CHRENT( I , J) 

11225  CONTINUE 

11230  CONTINUE 

DO  14344  1=1, LZMAX 

WTOT  =  WTOT  +  CHRW(I) 

WTABS  =  WTABS  +  CHRA(I) 

14344  CONTINUE 

VTESC  =  CHRI( LZMAX) 

WTCON  =  (WTABS  +  WTESC  +  WKILL)/WTOT 

TEMP(l)  =  WTCON 

CALL  STATS (TEMP , ISIG, 1 ,KPUT) 

WTCON  =  TEMP(l) 

C  CALCULATE  THE  AVERAGE  VALUES  OF  THE  TOTALS  AND  THEIR  ERRORS 
CALL  STATS ( CHRN , I S IG 1 , LZMAX , KPUT) 

CALL  STATS ( CHRW , ISIG2 , LZMAX , KPUT) 

CALL  STATS ( CHRE, I SIG3, LZMAX, KPUT) 

CALL  STATS ( CHRA , I S I G4 , LZMAX , KPUT ) 

CALL  STATS ( CHRO , ISIG5 , LZMAX , KPUT ) 

CALL  STATS (CHRI, I SIG6, LZMAX, KPUT) 

ENDIF 

C 

IF  ( IVIRS  .EQ.  1)  THEN 
DO  96  IR=1,NRVIR 
DO  86  IC=1 ,NCVIR 
DO  76  IA=1 .NAVIR 

SCOVOL( IA, IC , IR)=SCOVOL( IA , IC , IR) /CIMAX 
76  CONTINUE 

86  CONTINUE 

96  CONTINUE 

C 

ENDIF 

C 

*1 .OUTPUT. 1375 

IF( ICKV  .EQ.  1)  THEN 

C*****fr>v**iWrfr**7Wr>v*0}{ANGES  ADDED  HERE******************’*'"** 

NPUT=NKB I NS  - 1 
DO  11220  1=1, LZMAX 
WRITE(NPRT, 9000)  I 
9000  FORMAT( / , ' ZONE  =' ,13,/, 

1  ' WVLNTH ' ,T14, 'NUMBER' ,T22,'STAT' ,T33 ,' WEIGHT' ,T41,'STAT' 

2  T52, 'ENERGY' ,T60,'STAT' ,T7 1 ABSORB ' ,T79,'STAT' ,T90, 

3  '  LEAVE' ,T97, 'STAT' ,T108, '  ENTER' ,T115, ' STAT' ,/ , 


1 

2 

3 

11110 

11115 

11120 

11125 

11130 

11135 

11140 

11145 

11150 

11155 

11160 

11165 

11170 

11180 

1 

2 


’ . 1  ,T14, ' . 1  ,T22 , ' - 

T52  , ' . '  ,T60 ,T7 1 , '  - 

' . ’  ,T97,' - ’  ,T108,’ - 

DO  11110  J=1,NMBINS-1 
TEMP( J)=CHRNUM( I , J) 

CONTINUE 

CALL  STATS(TEMP, ISIGA, NPUT, KPUT) 
DO  11115  J=1 ,NMBINS- 1 
CHRNUM( I , J)=TEMP( J) 

CONTINUE 

DO  11120  J=1,NMBINS-1 
TEMP( J)=CHRWAT( I ,  J) 

CONTINUE 

CALL  STATS (TEMP, I SIGB , NPUT, KPUT) 
DO  11125  J=1,NMBINS-1 
CHRWAT( I , J)=TEMP( J) 

CONTINUE 

DO  11130  J=1 ,NMBINS-1 
TEMP( J)=CHRENG( I , J) 

CONTINUE 

CALL  STATS(TEMP , ISIGC ,NPUT,KPUT) 
DO  11135  J=1 ,NMBINS- 1 
CHRENG( I , J)=TEMP( J) 

CONTINUE 

DO  11140  J=1 ,NMBINS-1 
TEMP( J)=CHRABS( I ,  J) 

CONTINUE 

CALL  STATS ( TEMP , I S I GD , NPUT , KPUT) 
DO  11145  J=1,NMBINS-1 
CHRABS( I , J)=TEMP( J) 

CONTINUE 

DO  11150  J=1 ,NMBINS- 1 
TEMP( J)=CHROUT( I , J) 

CONTINUE 

CALL  STATS (TEMP, I SIGE, NPUT, KPUT) 
DO  11155  J=1 ,NMBINS-1 
CHROUT ( I , J ) =TEMP ( J ) 

CONTINUE 

DO  11160  J=1 ,NMBINS-1 
TEMP( J)=CHRENT( I , J) 

CONTINUE 

CALL  STATS ( TEMP , I S I GF , NPUT , KPUT ) 
DO  11165  J=1 .NMBINS-l 
CHRENT( I , J)=TEMP( J) 

CONTINUE 
DO  11170  J=1 , 17 

IWVL( J)=1600  +  J*200 
CONTINUE 

DO  11180  J=18 ,27 

1WVL( J)=5000  +  (J-17)*500 
CONTINUE 

DO  11190  J=1 .NMBINS- 1 


T33 , 1 --- 
,T79 
’ ,T1 15 , ’ 


,T41 , 1 - 

- - ’ ,T90, 


WRITE(NPRT,  11181)  IWVL(  J)  ,CHRNUM(  I ,  J) ,  ISIGV  J)  , 
CHRWAT( I , J) , 


ISIGB(  J) ,CHRENG( I , J) ,ISIGC( J) ,CHRABS( I , J) , ISIGD( J) , 


11190 

11220 

93867 


12181 


12190 


12191 


3  CHROUT(I,J),ISIGE(J),CHRENT(I,J),ISIGF(J) 

11181  F0RMAT(I5,T8,1PE12.4,T22,I3,T27,1PE12.4,T41,I3,T46, 

1  1PE12.4, 

2  T60 , 13 ,T65 , 1PE12.  4 ,T79 , 13 ,T84 , 1PE12. 4,T97 , 13 ,T102 , 

3  1PE12. 4.T115 ,13) 

11190  CONTINUE 

11220  CONTINUE 

WRITECNPRT, 93867) 

93867  FORMATC /,' TOTALS  FOR  ALL  ZONES',/, 

1  'ZONE' ,T14, 'NUMBER' ,T22, 'STAT* ,T33 , 'WEIGHT' ,T41, 'STAT' , 

2  T52, 'ENERGY' ,T60, 'STAT' ,T71 , 'ABSORB' ,T79, 'STAT' ,T90, 

3  '  LEAVE' ,T97,' STAT’ ,T108,'  ENTER* ,T115 ,' STAT' ,/ , 

1  ' . '  ,T14 , ' . '  ,T22 , ' - '  ,T33 , ' . '  ,141,' - 

2  T52 , ' . '  ,T60,' - '  ,T71 , ' . ’  ,T79,' - '  ,T90, 

3  ' . '  ,T97  , ' - '  ,T108 , ' . '  ,T115 , ' - '  ,/) 

DO  12190  J=1 ,LZMAX 

WRITECNPRT, 12181)  J,CHRN( J) ,ISIG1( J) , 

1  CHRW(J), 

2  ISIG2( J) , CHREC  J) ,ISIG3(J) ,CHRA( J),ISIG4(J), 

3  CHRO(J) ,ISIG5(J) ,CHRI( J) ,ISIG6(J) 

12181  FORMAT( I5,T8,1PE12. 4 ,T22 , 13 ,T27 , 1PE12.  4 ,T41 , 13 ,T46 , 

1  1PE12.4, 

2  T60 , 13 ,T65 , 1PE12.  4,T79 , 13 ,T84, 1PE12. 4 ,T97 , 13 ,T102 , 

3  1PE12. 4,T115 , 13) 

12190  CONTINUE 

WRITE CNPRT, 12191)  WTCON 

12191  F0RMAT(/,1X, 'THE  WEIGHT  CONSERVATION  FRACTION  IS  ',E12.4) 

END  IF 

C********?v************* END  OF  CHANGED  OUTPUT*********************** 
90100  FORMATC 13, 3X , 6 ( E 1 2 .  4,13)) 

IF( IVIRS  .EQ.  1)  THEN 
WRITECNPRT, 90200) 

90200  FORMAT( /' CERENKOV  VIRTUAL  SOURCE,  PHOTONS/SQCM/STR’ ) 

DO  6030  IVR=1 ,NRVIR 

WRITECNPRT, 90300)  RVIRB( IVR) ,RVIRB( IVR+1) 

90300  FORMATC'  R-MIN'  ,E12. 4/ '  R-MAX'  ,E12. 4,24X, 

1  '  UPPER  AND  LOWER  POLAR  COSINES  ') 

DO  6020  LVC=1 ,NCVIR , 5 
MVC=MIN0C  LVC+4 ,NCVIR) 

WRITECNPRT , 9040) C  CVIRBC IVC) , IVC=LVC ,MVC) 

WRITECNPRT, 9050)CCVIRBCIVC+1),IVC=LVC,MVC) 

9040  FORMATC 23X, 'MINIMUM' ,5X, 1  MAXIMUM’ ,6E15.4) 

9050  FORMATC 23X, ’AZIMUTH1 ,5X, ’AZIMUTH’ ,6E15.  4) 

DO  6010  IVA=1 ,NAVIR 
L=0 

DO  6000  IVC=LVC ,MVC 
L=L+1 

TEMP( L)=SCOVOLC IVA , IVC , IVR) 

6000  CONTINUE 

CALL  STATS(TEMP,ISIG,L,KPUT) 

WRITECNPRT, 9060)AVIRB( IVA) ,AVIRBC IVA+1) , 

1  ( TEMPC IVC) , ISIG( IVC) , IVC=1 ,L) 

9060  FORMAT ( 18X , 2E12. 4 , 6C E12. 2 , 13) ) 

6010  CONTINUE 

6020  CONTINUE 


90200 


90300 
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6030  CONTINUE 

C 

END  IF 

C 

*1, OUTPUT. 1383 

I F ( ICKV  .EQ.  1)  THEN 
WKILL  =  CZERO 
DO  112  J=1 . LZMAX 
DO  111  IJ=1,100 

CHRNUM( J, IJ)=CZERO 
CHRVAT( J, IJ)=CZERO 
CHRENG( J , I J)=CZERO 
CHRABS( J , I J)=CZERO 
CHROUT( J,IJ)=CZERO 
CHRENT( J , I J)=CZERO 
111  CONTINUE 

112  CONTINUE 
END  IF 

C 

I  F(  IV  IRS  .EQ.  1)  THEN- 
DO  115  IR=1 jNRVIR 
DO  114  IC=1 jNCVIR 
DO  113  IA=1 ,NAVIR 

SCOVOL( IA, IC , IR)=CZERO 

113  CONTINUE 

114  CONTINUE 

115  CONTINUE 

C 

END  IF 

C 

*1, INPUT. 1853 

9051  VRITE(*,9070) 

9070  FORMAT 

1  (  GFRF.MATURE  EOF  WHILE  READING  CKV  OR  VIRS  DATA') 

GO  TO  99.-9 

9052  VRITE( *,9071) 

9071  FORMAT 

1  ( ' OILLEGAL  DATA  IN  FIELD  WHILE  READING  CKV  OR  VIRS  DATA') 

GO  TO  9999 

C 

*1  .EHIST. 522 

IF  (ICKV  .EQ.  1)  THEN 

C  ***  SAVE  ZONE  INDEX  FOR  POSSIBLE  CALL  TO  TRANSRAD 
IRSAV  =  IR 
IRPSAV  =  IRPRIM 
DISADD=DISTO 

I F( MARK  . NE.  1)  DISADD=DIST  +  l.D-07 
XCHR=X  +  DISADD*STH( 1)*CPH( 1) 

YCHR=Y  +  DISADD*STH( 1)-SPH( 1) 

ZCHR=Z  +  DISADD“CTH( 1 ) 

CALL  CKVRAD 

1  ( XCHR , YCHR , ZCHR , CTH , STH , CPH , SPH , W ,T) 

C  ***  REINSTATE  ZONE  INDEX  FOR  POSSIBLE  CALL  TO  TRANSRAD 
IR  =  IRSAV 
IRPRIM  =  IRPSAV 


IF( IRPRIM  .  LT.  0)  RETURN 


END  IF 
*I,DISTA.  133 

NBOSAV  =  MA(LSAVE) 

*D, PAREM. 6 

$  KLOOP,  LOOP,  ITYPE ,  NBOSAV 
*1, KNOCK.  367 

C  *CERENKOV  PRODUCTION  SELECTED  RANDOMLY  ALONG  A  SUBSTEP* 

C  ******************************************************* 

SUBROUTINE  CKVRAD 

1  ( XS , YS , ZS , CPOL, SPOLjCAZI , SAZI , WBETA ,EBETA) 


VARIABLE  DESCRIPTIONS 


ACCEPT 

ALPHA 
ANGTCM 
C(  1) 

C(  2) 

C(  3) 
CAFT 
CAZI 


FRACTION  OF  CERENKOV  RADIATION  ACCEPTED  INTO 
ZONE 

1./ 137.0  (FINE  STRUCTURE  CONSTANT) 
l.E-8  (AN  ANGSTROM  IN  CENTIMETERS) 

X  PRESENT  DIRECTION  COSINE 

Y  PRESENT  DIRECTION  COSINE 

Z  PRESENT  DIRECTION  COSINE 

COSINE  AFTER  REFLECTION 

COSINE  OF  AZIMUTHAL  SCATTERING  ANGLE 


C**********CERENKOV  INPUTS******* 


CHER(l.J) 


CHER( 2 , J) 
CHERf  3 , J) 

CHER(4,J) 

CHER( 5 , J) 


MAl'ERIAL  NUMBER  FOR  ZONE  J 

1  =  VACUUM 

2  =  AIR 

3  =  C02 

4  =  GLASS 

5  =  QUARTZ 

EXTINCTION  COEFFICIENT  (1/CM) 

FRACTION  OF  LIGHT  TO  BE  MIRROR  REFLECTED 
FROM  ZONE  WALLS 

FRACTION  OF  LIGHT  TO  BE  DIFFUSE  REFLECTED 
FROM  ZONE  WALLS 

FRACTION  OF  LIGHT  TO  BE  ACCEPTED  FOR 
ENTRANCE  INTO  ZONE 


C Vr**** CERENKOV  OUTPUTS*** 
C 

C  CHRNUM( J,K)  NU 


CHRWAT( J,K) 


CHRENG( J,K) 


CHRABS(J.K) 


CHROUT(J.K) 


NUMBER  OF  CERENKOV  PHOTONS  PRODUCED  IN  ZONE  J 
FROM  WAVELENGTH  BIN  K 
PER  INCIDENT  SOURCE  PARTICLE 
CERENKOV  WEIGHT  PRODUCED  IN  ZONE  J 
FROM  WAVELENGTH  BIN  K  PER  INCIDENT 
SOURCE  PARTICE  WEIGHT  OF  1.0 
CERENKOV  ENERGY  PRODUCED  IN  ZONE  J 
FROM  WAVELENGTH  BIN  K  PER  INCIDENT 
SOURCE  PARTICLE 
=  CHRENG(I)  +  EP*WP 

CERENKOV  WEIGHT  ABSORBED  (KILLED)  IN  ZONE  J 
FROM  WAVELENGTH  BIN  K  PER 
INCIDENT  SOURCE  PARTICLE 
=  CHRABS(J')  +  WNEW  OR  WABS 
CERENKOV  WEIGHT  THAT  EXITS  ZONE  J 
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o  o  o  n  o  o  n 


FROM  WAVELENGTH  BIN  K  PER 
INCIDENT  SOURCE  PARTICLE 
=  CHROUT(J)  +  WP 

CHRENT( J,K)  CERENKOV  WEIGHT  THAT  ENTERS  ZONE 

FROM  WAVELENGTH  BIN  K  PER 
INCIDENT  SOURCE  PARTICLE 
=  CHRENT( J+l)  +  WP 


(J+l) 


C  CN( J) 

C  CP(J) 

C  CP2 
C  POL 
C  CPRE 
C  CPRO 
C  CSQ 
C 

C  CX 
C  CZ 
C  DIFFU 
C  EBETA 
C  EE 
C  EP 
C  FST 
C 

n 

c  J 

C  KALTOT 
C  KLIDE 
C  NR 
C  NRC 
C  NR  PC 
C  NRMAX 
C  NSC 
C  NZON 
C  PI 
C  PLANKC 
C  RA 

C  RANF(ARG) 

C 

C  RVIR 
C  SAFT 
C  SAZI 

C  SCOVOL( J, J, J) 
C  SHD 
C  SIG 
C  SPECU 
C  SPOL 
C  ST 
C  oTP 
C  TIME 
C  V 
C  V2 
C  V2N2 

C  VARIND( I , J) 

C 

C  VOLVIR 


NORMAL  FROM  THIS  REGION  TOWARDS  NEXT  REGION 
DIRECTION  COSINES  IN  MOVING  FRAME 
COSINE  SQUARED 

COSINE  OF  POLAR  SCATTERING  ANGLE 
SUM  OF  SQUARES  OF  CURRENT  DIRECTION  COSINES 
PRODUCTION  COSINE  IN  THE  MOVING  FRAME 
SQUARE  ROOT  OF  SUM  OF  SQUARES  OF  X  AND  Y  DIRECTION 
COSINES 

X  DIRECTION  COSINE 
Z  DIRECTION  COSINE 
PROBABILITY  FOR  DIFFUSE  REFLECTION 
ELECTRON  ENERGY 
GAMMA  SQUARED 
ENERGY  OF  CERENKOV  PHOTON 
FRACTION  OF  LIGHT  TO  BE  TRANSMITTED  OUT 
THROUGH  ZONE  BOUNDARY 
=  1.0  -  CHER( 3 , J)  -  CHER(4 , J) 

COUNTER 

MAXIMUM  NUMBER  OF  ALLOWED  REFLECTIONS 

COUNTER  FOR  NUMBER  OF  REFLECTIONS 

PRESENT  USER  REGION 

PRESENT  CODE  REGION 

NEXT  CODE  REGION 

MAX  ZONE  NUMBER 

SURFACE  CROSSED  BY  RAY 

ZONE  NUMBER 

3.  14159 

1. 2405E- 10 

RANDOM  NUMBER 

RANDOM  NUMBER  PRODUCING  FUNCTION 

ARG  =  0  PRODUCES  UNIFORM  RANDOM  NUMBER 
VIRTUAL  SOURCE  RADIUS 

SINE  OF  AZIMUTHAL  SCATTERING  ANGLE 

DATA  FOR  VIRTUAL  SOURCE 

LENGTH  OF  ELECTRON  SUBSTEP  IN  CM 

PRODUCTION  OF  CERENKOV  PHOTONS  PER  CM  PER  MEV 

PROBABILITY  FOR  SPECULAR  REFLECTION 

SINE  OF  POLAR  SCATTERING  ANGLE 

DISTANCE  TO  BOUNDARY 

SAMPLED  LENGTH  ALONG  ELECTRON  SUBSTEP  (CM) 

RATIO  OF  ELECTON  VELOCITY  TO  SPEED  OF  LIGHT  (BETA) 
BETA  SQUARED 

BETA  SQUARED  TIMES  INDEX  OF  REFRACTION  SQUARED 
INDEX  OF  REFRACTION  FOR  MATERIAL  I  AND  WAVELENGTH 
BIN  J 

VOLUME  FOR  VIRTUAL  SOURCE 
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C  WABS 
C  WVLNTH(I) 
C 
C 
C 

C  WBETA 
C  WNEW 
C  WP 
C  WT(J) 

C  WZERO 
C  XB( J) 

C  XCHR(J) 

C  XPI(J) 

C  XPO(J) 

C  XREF 
C  XS 
C  XVIR 
C  YREF 
C  YS 
C  YVIR 
C  2E 
C  ZS 
C 


CERENKOV  WEIGHT  ABSORBED  IN  ZONE 
WAVELENGTH  BIN  I 

THESE  BINS  RUN  FROM  1800  TO  5000  ANGSTROMS 
IN  200  ANGSTROM  WIDE  BINS  AND  THEN  FROM 
5500  TO  10000  ANGSTROMS  IN  500  ANGSTROM  BINS 
ELECTRON  WEIGHT 

REDUCED  WEIGHT  OF  CERENKOV  PHOTON 
WEIGHT  OF  CERENKOV  PHOTON 
TEMPORARY  POSITION  ARRAY 
J.NTTIAL  WEIGH.  OF  CERENKOV  PHOTON 
TEMPORARY  POSITION  ARRAY 
PRESENT  POSITION 

POSITION  AT  BOUNDARY  JUST  INSIDE  THIS  REGION 
POSITION  AT  BOUNDARY  JUST  INSIDE  NEXT  REGION 
INDEX  OF  REFRACTION 

CURRENT  X  POSITION  (INPUT  IN  PARAMETER  LIST) 

VIRTUAL  SOURCE  X  POSITION 

INDEX  OF  REFRACTION  FOR  NEXT  ZONE 

CURRENT  Y  POSITION  (INPUT  IN  PARAMETER  LIST) 

VIRTUAL  SOURCE  Y  POSITION 

CHARGE  ON  SOURCE  PARTICLE  IN  ELECTRON  VOLTS 
CURRENT  Z  POSITION  (INPUT 
IN  PARAMETER  LIST) 


C********************************************************************* 


C 

*CA  CNSTNT 
*CA  PARAMS 
*CA  OUT 


*CA  CALC 
*CA  PAREM 
*CA  STTS 
*CA  GOMLOC 
*CA  DBG 
*CA  ORGI 
*CA  COG 
*CA  SECZN 

DOUBLE  PRECISION  DSEED,  IRAN 
COMMON  /VAXRAN/  IRAN 

DIMENSION  C( 3) ,CN( 3 ) ,CP( 3) ,XCHR( 3) ,XPI( 3) ,XPO( 3) 
DIMENSION  CDF(20, 100)  ,PDF( 100) 

DATA  PLANKC , ANGTCM/ 1 .  2405D-10, 1.  D-8/ 

C  RANDOM  NUMBER  STATEMENT  FUNCTION  FOR  IBM 
RANF( DSEED )=  GGUBFS( DSEED) 

C  SET  CURRENT  NUMBER  OF  MATERIALS  WITH  DATA  FOR  INDEX 
C  OF  REFRACTION  IN  TABLE  ABOVE. 


NMTRLS  =  5 

C  COUNT  THE  NUMBER  OF  TIMES  THIS  ROUTINE  IS  CALLED 
NCALL  =  NCALL  +  1 
C  SET  SOME  CONSTANTS 

ALPHA  =  1. DO/137.  DO 
PI=DACOS( - 1. DO) 

TIME=1. DO 


M=1 


C  CALCULATE  DIRECTION  COSINES  FOR  SOURCE  PARTICLE 


C(  1 )  =  SP0L'vCAZI 
C( 2)  =  SPOL*SAZI 
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C( 3)  =  CPOL 

C  SAMPLE  RANDOM  LOCATION  FOR  CERENKOV  PRODUCTION  ALONG  A  SUBSTEP 
STP  =  SHD*RANF ( I RAN ) 

C  SET  ELECTON  LOCATION 
XCHR(l)  =  XS 
XCHR( 2)  =  YS 
XCHR( 3)  =  ZS 

C  PRESENT  ZONE  OF  ELECTRON 
NRC  =  LBCZ 
NR  =  IROR(NRC) 

NRMAX  =  NZON 

C  STEP  BACK  PART  WAY  ALONG  SUBSTEP  TO  SAMPLE  CERENKOV  PRODUCTION 
DO  100  J=1 , 3 
VT(J)  =  C(J) 

XB( J)  =  XCHR(J)  +  C( J)*( STP  -  SHD) 

XCHR(J)  =  XB( J) 

100  CONTINUE 

C  FIND  ZONE  OF  CERENKOV  PRODUCTION  AND  COMPUTE  UNCOLLIDED 
C  ESCAPE  DISTANCE. 

CALL  ZONEA 

NEXT  ZONE  TO  BE  ENCOUNTERED  BY  CERENKOV  RAY. 

NRC=IRPRIM 
C  PRESENT  ZONE  OF  CERENKOV  PRODUCTION 
NR=IR 

C  WHAT  MATERIAL  IS  ELECTRON  IN  FOR  CERENKOV  PRODUCTION? 

MTRL= I D I NT( CHER ( 1 , NR ) ) 

C  GET  THE  INDEX  OF  REFRACTION  FOR  THAT  MATERIAL 
XREF=VARIND(MTRL, 27) 

C  ***  CALCULATE  VELOCITY  OF  SOURCE  PARTICLE  *** 

C  ELECTRON  ENERGY  NORMALIZED  TO  REST  MASS  ENERGY 
C  TOTAL  ENERGY  =  KINETIC  ENERGY  +  REST  MASS  ENERGY  (EE  =  GAMA) 

EE  =  1. DO  +  EBETA/. 511D0 
C  NORMALIZED  VELOCITY  (V2  =  BETA  SQUARED) 

V2  =  l.DO  -  1. DO/(EE*EE) 

C  (V  =  BETA) 

V  =  DSQRT(V2) 

C  (BETA  SQUARED  *  INDEX  OF  REFRACTION  SQUARED) 

V2N2  =  V2*XREF*XREF 

C  CALCULATE  THE  CDF  FOR  ALL  THE  MATERIALS  AVAILABLE 
IF(NCALL  .EQ.  1)  THEN 

CALL  CKVCDF( V , NMTRLS , VARIND , WVLNTH , CDF) 

END  IF 

C  ***  DOES  PARTICLE  MEET  CERENKOV  THRESHOLD  *** 

IF( XREF  .LE.  l.DO)  GO  TO  1600 
IF( V2N2  .LE.  l.DO)  GO  TO  1600 

C  THE  CERENKOV  PHOTON  WEIGHT  TO  BE  USED  IS  THE  TOTAL  INTEGRATED 
C  WEIGHT  OVER  THE  ENTIRE  RANGE  OF  WAVELENGTHS  AND  FOR  THE 
C  ENTIRE  LENGTH  OF  THE  SUBSTEP 
WP  =  0.  DO 
DO  150  1=2 ,NMBINS 

WP  =  2.  DO*PI*ALPHA*SHD*WBETA*DABS(( 1. DO/(ANGTCM*WVLNTH( 1-1) ) 

1  -  1 .  DO / ( ANGTC M*WV LNTH ( I ) ) ) * 

2  (1.  DO-1.  DO/V/V/VARIND(MTRL, I- 1)/VARIND(MTRL, I - 1 ) ) ) 

3  +  WP 
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150  CONTINUE 

C  SAMPLE  TO  GET  CORRECT  CERENKOV  BIN 
CALL  CKWVL(CDF ,MTRL,NCBIN) 

C  ***  GET  INDEX  OF  REFRACTION  *** 

XREF  =  VARIND(MTRL,NCBIN) 

C  ***  is  INDEX  OF  REFRACTION  ABOVE  CUTOFF  FOR  CERENKOV?  *** 

IF(XREF  .  LE.  1. DO)  GO  TO  1600 
C  CALCULATE  PHOTON  ENERGY 

EP  =  PLANKC/(ANGTCM*WVLNTH(NCBIN)) 

C  ***  CALCULATE  COSINE  FOR  CERENKOV  PHOTON  IN  MOVING  FRAME 
C  ***  =  1  /  (BETA  *  INDEX  OF  REFRACTION) 

C  ***  =  RATIO  OF  SPEED  OF  LIGHT  IN  MEDIUM  TO  ELECTRON  VELOCITY 

CPRO  =  1 . DO/ ( V*XREF ) 

C  ***  TALLY  PHOTON  NUMBER  COUNT,  WEIGHT,  AND  ENERGY  PRODUCED  HERE 
CHRNUM(NR,NCBIN)  =  CHRNUM (NR,NCBIN)  +  1.D0 
CHRWAT(NR.NCBIN)  =  CHRWAT(NR,NCBIN)  +  WP 
CHRENG(NR.NCBIN)  =  CHRENG(NR.NCBIN)  +  EP*WP 
C  ***  BRING  POLAR  ANGLE  BACK  INTO  LAB  FRAME  *** 

CALL  FOLD( CPOL , SPOL, CAZ I , SAZI , CPRO , CZ , SZ , CX , SX) 

C  ***  USE  SAME  INITIAL  PHOTON  WEIGHT 
WZERO  =  WP 

c  ***  SET  DIRECTION  COSINES  AND  PREPARE  FOR  RAY  TRACE  *** 

C(l)  =  SZ*CX 
C( 2)  =  SZ*SX 
C( 3 )  =  CZ 

C  ***  CONDUCT  RAY  TRACE  OF  CERENKOV  PHOTON  THROUGH  GEOMETRY  *** 

NUMCT  =  NUMCT  +  1 
C  SET  INITIAL  COLLISTION  COUNTER 
KLIDE  =  0 

C  ***  NR,  PRESENT  (USER)  REGION 
C  ***  NRC,  PRESENT  (CODE)  REGION 
C  *•'*  NRPC ,  NEXT  (CODE)  REGION 
C  ***  X,  PRESENT  POSITION  OF  CERENKOV  RAY 
C  ***  C,  PRESENT  DIRECTION  OF  CERENKOV  RAY 
C  ***  ST,  DISTANCE  TO  BOUNDARY 
C  ***  NRP,  NEXT  (USER)  REGION 
C  ***  NSC,  SURFACE  ENCOUNTERED  BY  RAY 

C  ***  XPI,  POSITION  OF  RAY  AT  BOUNDARY  JUST  INSIDE  THIS  REGION 

C  ***  XPI,  POSITION  OF  RAY  AT  BOUNDARY  JUST  INSIDE  THE  NEXT  REGION 

C  ***  CN,  NORMAL  DIRECTION  FROM  THIS  REGION  TOWARD  NEXT  REGION 

700  CALL  PATH( NR , NRC ,XCHR ,C , ST, NRPC , NSC , XPI ,XPO ,CN) 

C  ***  CHECK  FORMATTING  HERE  *** 

C  ***  IF  NO  ZONE  IS  FOUND  THEN  PHOTON  ESCAPES 
IF( IRPRIM  .  LT.  0)  THEN 
M  =  5 

GO  TO  1600 
END  IF 
M  =  3 

C  ***  COMPUTE  THE  PHOTON  WEIGHT  ATTENUATION  AND  DETERMINE  THE  WEIGHT 
C  ***  LEFT  FOR  FURTHER  RAY  TRACE. 

WNEW  =  WP*DEXP( -ST*CHER( 2 ,NR) ) 

C  ***  UPDATE  THE  WEIGHT  (CALCULATE  THE  WEIGHT  ABSORBED  IN  ZONE) 

WABS  =  WP  -  WNEW 
WP  =  WNEW 

C  ***  TALLY  THE  ABSORBED  WEIGHT  FOR  THAT  ZONE 


CHRABS(NR,NCBIN)  =  CHRABS(NR.NCBIN)  +  WABS 

C  ***  THIS  SECTION  CREATES  TALLIES  FOR  VIRTUAL  SOURCE 
IF( IVIRS  .EQ.  1)  THEN 

ICVIR  =  IBFIND(C( 3) ,NCVIR,CVIRB) 

IF( ICVIR  .LE.  0)  GO  TO  1000 
SVIR  =  (HVIR  -  XCHR( 3) )/C( 3) 

XVIR  =  XCHR(l)  +  SVIR*C(1) 

YVIR  =  XCHR( 2)  +  SVIR*C( 2) 

RVIR  =  DSQRT(XVIR*XVIR  +  YVIR*YVIR) 

IRVIR  =  IBFIND(RVIR,NRVIR ,RVIRB) 

IF( IRVIR  .LE.  0)  GO  TO  1000 
CTVIR  =  1. DO 
STVIR  =  0.  DO 
IFCRVIR  .  GT.  0. DO)  THEN 
CTVIR  =  XVIR/RVIR 
STVIR  =  YVIR/RVIR 
END  IF 

AVIR  =  0. DO 

CSQ  =  DSQRT(C(1)*C(1)  +  C(2)*C(2)) 

IF(  CSQ  .LE.  0.  DO)  GO  TO  900 

AVIR  =  (CTVIR*C(1)  +  STV I R*C ( 2 ) ) /CSQ 

IF(DABS(AVIR)  . GT.  l.DO)  AVIR  =  AVIR/DABS( AVIR) 

AVIR  =  DACOS(AVIR) 

900  IAVIR  =  IBFIND( AVIR ,NAVIR , AVIRB) 

IF( IAVIR  .LE.  0)  GO  TO  1000 

VOLVIR  =  2.  D0*PI*(RVIRB( IRVIR  +  1)**2-RVIRB( IRVIR)**2) 

1  *(CVIRB( ICVIR  +  1)  -  CVIRB( ICVIR)) 

2  *(AVIRB( IAVIR  +  1)  -  AVIRB( IAVIR) ) 

SCOVOL( IAVIR , ICVIR , IRVIR)=SCOVOL( IAVIR, ICVIR , IRVIR)+WP/VOLVIR 
END  IF 

C************  END  OF  VIRTUAL  SOURCE  SECTION  ***************** 

C  ***  INCREMENT  THE  NUMBER  OF  REFLECTIONS  (COLLISIONS  OF  PHOTON) 

1000  KLIDE  =  KLIDE  +  1 

C  ***  PUT  AN  UPPER  LIMIT  ON  THE  NUMBER  OF  REFLECTIONS  ALLOWED 

C  ***  (IE.  DON'T  GET  IN  AN  INFINITE  LOOP  BETWEEN  WALLS) 

IF(  (KLIDE  .GT.  KALTOT)  .OR.  (WP  .  LT.  WZER0*1.  D-3) )  THEN 
WKILL  =  WKILL  +  WP 
GO  TO  1600 
END  IF 

C  ***  INITIALIZE  SPECULAR  AND  DIFFUSE  REFLECTION  COEFFICIENTS 
SPECU  =  0. DO 
DIFFU  =  0. DO 

C  ***  SET  NEXT  REGION  THAT  RAY  WILL  PASS  INTO 
NRP  =  NRPC 

C  ***  CHECK  TO  SEE  IF  WE  HAVE  A  NEXT  REGION 

c  ***  if  SO  THEN  GET  THE  TRUE  SPECULAR  AND  DIFFUSE  COEFFICIENTS 
IF ( NRPC  .GT.  0)  NRP=IROR( NRPC) 

IF(NRP  .GT.  0)  THEN 
SPECU  =  CHER( 3 ,NRP) 

DIFFU  =  CHER (4, NRP) 

END  IF 

C  ***  CHECK  TO  SEE  IF  PHOTON  UNDERGOES  SPECULAR  OR  DIFFUSE  REFLECTION 
RA  =  RANF( IRAN) -SPECU 
IF(  RA  .LE.  0.  DO)  THEN 

C  ***  SPECULAR  REFLECTION  *** 


SI 


CPRE  =  C(1)*CN(1)  +  C(2)*CN(2)  +  C(3)*CN(3) 

1100  PSI  =  2.  D0*CPRE 

C  ***  NEW  DIRECTION  =  OLD  DIRECTION  -  2*( NORMAL  COMPONENT) 

C  ***  VECTR2  =  VECTR1  -  2*(NMLVECT  DOT  VECTR 1 ) *NMLVECT 
C( 1)  =  C( 1)  -  PSI*CN(1) 

C(2)  =  C C 2 )  -  PSI*CN(2) 

C( 3)  =  C( 3)  -  PSI*CN( 3) 

C  ***  SET  NEW  POSITION  OF  THE  RAY 
DO  1200  J=1 , 3 

XCHR(J)  =  XPI(J) 

1200  CONTINUE 

C  ***  GO  BACK  TO  PATH  AND  FOLLOW  THE  RAY  THROUGH  THIS  NEXT  ZONE 
GO  TO  700 
END  IF 

C  ***  CHECK  FOR  DIFFUSE  REFLECTION  *** 

RA  =  RA  -  DIFFU 
IF(RA  .LE.  0.  DO)  THEN 

C  ***  SAMPLE  RANDOM  POLAR  COSINE  RELATIVE  TO  SURFACE  NORMAL 
CP( 3)  =  -RANF( IRAN) 

C  ***  SINE  OF  POLAR  SCATTERING  ANGLE 

SINE  =  DSQRT( 1.  DO  -  CP(3)*CP(3)) 

C  ***  AZIMUTHAL  SCATTERING  ANGLE  RANDOM  BETWEEN  +  OR  -  PI 
THE  =  PI*( 2.  D0*RANF( IRAN)  -  1.D0) 

CP(1)  =  SINE*DCOS(THE) 

CP( 2 )  =  SINE*DSIN(THE) 

C  ***  ROTATE  DIRECTION  COSINES  BACK  INTO  LAB  FRAME 
CALL  ROTCPC(CN,CP,C) 

C  ***  UPDATE  CURRENT  POSITION  OF  THE  RAY 
DO  1300  J=1 ,3 

XCHR(J)  =  XPI(J) 

1300  CONTINUE 

C  ***  GO  BACK  TO  PATH  AND  RAY  TRACE  THROUGH  THE  NEXT  ZONE 
GO  TO  700 
END  IF 

C  ***  IF  THEFT,  TS  NO  REFLECTION,  TRANSMIT  THE  CERENKOV  PHOTON  *** 

C  ***  SET  COSINE  OF  INCOMING  RAY  RELATIVE  TO  SURFACE  NORMAL 
CPRE  =  C(1)*CN(1)  +  C( 2)*CN( 2)  +C(3)*CN(3) 

YREF  =  1. DO 

C  ***  SET  ZONE  NUMBER  TO  THAT  OF  THE  ZONE  THE  RAY  IS  ENTERING 
NRP  =  NRPC 

IF( NRPC  . GT.  0)  NRP=IROR(NRPC) 

ACCEPT  =  1.  DO 

C  ***  ip  THERE  IS  A  NEXT  ZONE  SET  INDEX  OF  REFRACTION  AND  FRACTION 
C  ***  ACCEPTED 

IF(NRP  .GT.  0)  THEN 

MTRL1  =  IDINT(CHER( 1 ,NRP) ) 

YREF  =  VARIND(MTRI.I.NCBIN) 

ACCEPT  =  CHER( 5 , NRP ) 

END  IF 

C  ***  TEST  TO  SEE  IF  THE  PHOTON  IS  REFLECTED  OR  TRANSMITTED 
IF(RANF( IRAN)  .  GT.  ACCEPT)  GO  TO  1100 
C  ***  DOES  SNELLS  LAW  ALLOW  TRANSMITION?  (SNELL'S  LAW  IN  TERMS  OF 
C  ***  COSINE  OF  ANGLE  SQUARED) 

CP2  =  l.DO  -  XREF*XREF*( 1 . DO  -  CPRE*CPRE) /( YREF* YREF) 

C  ***  IF  COSINE  OF  ANGLE  SQUARED  LESS  THAN  OR  EQUAL  ZERO  THEN 
C  ***  SPECULAR  REFLECTION  OCCURS  (OUTSIDE  CRITICAL  ANGLE) 


# 


82 


IF( CP2  .  LE.  0. DO)  GO  TO  1100 
SPRE  =  DSQRT( 1.  DO  -  CPRE*CPRE) 

CAFT  =  DSQRT(CP2) 

SAFT  =  DSQRT( 1.  DO  -  CP2) 

SRAT  =  SAFT/ SPRE 
CRAT  =  CAFT  -  SRAT*CPRE 

C  ***  DIRECTION  COSINES  OF  RAY  AFTER  ENTERING  NEW  ZONE  IN  LAB  FRAME 
DO  1400  J=1 ,3 

C(J)  =  SRAT*C(J)  +  CRAT*CN(J) 

1400  CONTINUE 

C  ***  INCREMENT  CERENKOV  WEIGHT  EXITING  PREVIOUS  ZONE 
CHROUT(NR,NCBIN)  =  CHROUT(NR.NCBIN)  +  WP 
C  ***  IS  THERE  A  NEXT  ZONE? 

IF( NRP  .GE.  0)  THEN 

C  ***  INCREMENT  CERENKOV  WEIGHT  ENTERING  NEXT  ZONE 
CHRENT( NRP ,NCBIN)  =  CHRENT( NRP ,NCBIN)  +  WP 
C  ***  UPDATE  ZONE  AND  INDEX  OF  REFRACTION 
NR  =  NRP 
N’RC  =  NRPC 
XREF  =  YREF 

C  ***  SET  RAY  POSITION  JUST  INSIDE  NEXT  ZONE  BOUNDARY 
DO  1500  J=1 , 3 

XCHR(J)  =  XPO(J) 

1500  CONTINUE 

C  ***  GO  BACK  TO  PATH  AND  FOLLOW  THE  RAY  THROUGH  THIS  NEXT  ZONE 
I F ( NRP  .LT.  NRMAX )  GO  TO  700 
END  IF 

1600  RETURN 
END 

C  ***  SUBROUTINE  TO  ROTATE  COORDINATES  BACK  INTO  LAB  FRAME 
SUBROUTINE  ROTCPC(CP,D,C) 

IMPLICIT  REAL* 8  (A-H.O-Z) 

DIMENSION  CP( 3) ,D( 3) ,C( 3) ,E3(3),E2(3) ,E1(3) 

DO  100  J=1 , 3 

E  3 ( J )  =  CP(J) 

100  CONTINUE 

E22  =  E3(1)*E3(1)  +  E3(2)*E3(2) 

IF(E22  .LE.  0.  DO)  THEN 
E2f 1)  =  1. DO 
E2( 2)  =  0. DO 
ELSE 

E22  =  DSQRT(E22) 

E2( 1 )  =  -E3(2)/E22 
E2( 2)  =  E3( 1)/E22 
END  IF 

E 2 (  3 )  =  0.  DO 

El(’)  =  E2( 2)*E3( 3 )  -  E2(3)*E3(2) 

E 1( 2)  =  E2(3)*E3(1)  -  E2(1)*E3(3) 

El( 3)  =  E2(1)*E3(2)  -  E2(2)*E3(1) 

DO  200  J=1 , 3 

C(J)  =  D(1)*E1(J)  +  D(2)*E2(J)  +  D(3)*E3(J) 

200  CONTINUE 
RETURN 
END 

SUBROUTINE  PATHf  NR , NRG , X , C  , ST , N’RGP , NSCP  , XP I , XPO  , CN ) 

C  ***  RAY  TRACING  SUBROUTINE  **** 
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*CA  CNSTNT 
*CA  PARAMS 
*CA  PAREM 
*CA  STTS 
*CA  GOMLOC 
*CA  DBG 
*CA  ORGI 
*CA  SECZN 
*CA  COG 

COMMON/NRML/UN(  3 ) 

DIMENSION  X( 3) ,C(3) ,XPO(3) ,XPI(3) ,CN(3) 

C  ***  ACCEPT  VERSION 
DIST  =  CZERO 
DISTO  =  PINT 
IRPRIM  =  NRG 
IR  =  NR 

KLOOP  =  KLOOP  +  1 
DO  100  1=1,3 
XB( I )=X( I ) 

VT ( I ) =C ( I ) 

100  CONTINUE 
ICALL  =  0 
IZE  =  0 
CALL  DISTA(S) 

IF( IRPRIM  .GE.  0)  THEN 
ST  =  DIST 
DO  200  1=1,3 

XPI(I)  =  X(I)  +  (ST*0. 99999D0)*C(I) 

XPO(I)  =  X(I)  +  (ST*1.  0000 IDO )*C( I) 

200  CONTINUE 

NRGP  =  IRPRIM 

NSCP  =  LSURF 

CALL  NORMLC  NRG , LSAVE ) 

C  ***  SET  NORMAL  TO  DIRECTION  UNTIL  NORMAL  IS  AVAILABLE 
DO  300  1=1,3 

CN( I )  =  -UN( I ) 

300  CONTINUE 

END  IF 
RETURN 
END 

SUBROUTINE  NORML( IRR.NASC) 

C  ***  USER  BEWARE:  MUCH  OF  THIS  ROUTINE  IS  UNTESTED. 

C  ***  THAT  XB+WT*DIST  IS  ON  SURFACE  LSURF  OF  BODY  NASC. 

C  ***  CALL  MUST  BE  PROCEEDED  BY  A  CALL  TO  G1  OR  LOOKZ. 

*CA  CNSTNT 

*CA  PARAMS 

*CA  PAREM 

*CA  GOMLOC 

*CA  DBG 

*CA  ORGI 

*CA  SECZN 

*CA  COG 

COMMON  /NRML/UN( 3) 

DIMENS I ON  XP( 3 ) , VP( 3 ) , X( 3 ) , H( 3 ) , IB0X1 ( 3 ) , IB0X2 ( 3 ) 
DATA  IBOXl/10,7 ,4/ 

DATA  IB0X2/4 ,10,7/ 


IT  IS  ASSUMED 
THUS  A  NORMAL 
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LSUR  =  IABS(LSURF) 

DO  100  1=1,3 

XP(I)  =  XBCI)  +  DIST*WT(I) 

100  CONTINUE 

N  =  LOCRFG(IRR) 

NUM  =  NUMBOD( IRR)*4  +  N  -  4 
DO  200  I=N,NUM,4 

IF(NASC  .  EQ.  I)  GO  TO  400 
200  CONTINUE 

WRITE( 6 , 300)  IRR ,NASC , ( MA( I ) , I=N,NUM,4) 

300  FORMAT( IX, ' INVALID  REGION  OR  BODY  IN  NORMAL.  IR  =  ,17,  NA SC 
1  17,/, (1017)) 

CALL  PR( 1) 

STOP 

400  L  =  MA( 1+1) 

I TYPE  =  MA(L+1) 

K  =  MA( L+4) 

C  ARB  SPH  RCC  REC  TRC  ELL  BOX  WED  RPP 

GO  TO  (500,1200,1400,1500,1400,1600,1700,1500,1900),  I TYPE 
C  ***  ARB 

500  I=K  +  1  +  4*( LSUR  -  1) 

DO  1110  J=l,3 

UN(J)  =  FPD( I+J) 

1110  CONTINUE 

GO  TO  2000 
C  ***  SPH 

1200  DO  1210  1=1,3 

UN( I )  =  (XP(I)  -  FPD(K+I+l))/FPD(X+5) 

1210  CONTINUE 

GO  TO  2000 
C  ***  TRC  AND  RCC 
1400  H2=0.  DO 

DO  1410  1=1,3 

H( I )=FPD(K+4+I ) 

H2=H2  +  H( I)*H( I) 

1410  CONTINUE 

GO  TO  (1415,1415,1435),  LSUR 
1415  H2  =  DSQRT(H2) 

DO  1420  1=1,3 

UN ( I )  =  H(I)/K2 
1420  CONTINUE 

GO  TO  2000 
1435  DO  1440  1=1,3 

X(I)  =  XP(I)  ~  FPD(X+1+I) 

1440  CONTINUE 

XMU  =  0.  DO 
A2  =  0. DO 
DO  1445  1=1,3 

J  =  MOD( 1,3)  +  1 

A2  =  A2  +  (H(I)*X(J)  -  h(J)*X(I))**2 
XMU  =  XMU  +  FPD(K+4+I )*X( I ) 

1445  CONTINUE 

R  =  FPD(K+8)  -  FPD( K+9 ) 

IF( ITYPE  .EQ.  3)  R  =  0.  DO 
B 1  =  R/DSQRT(H2*(H2  +  R*R)) 

B2  =  l.D0/DSQRT(A2*(H2  +  R*R) ) 


S5 


1450 

1500 

1501 

C  *** 
1600 


1620 

1630 

1635 

1650 

1651 


C  *** 
1700 

1720 


1730 

1735 

C  *** 
1900 


DO  1450  1=1,3 

UN( I)  =  B2*(X( I)*H2  -  I1(I)*XMU)  +  B1*H(I) 

CONTINUE 

GO  TO  2000 

WRITE(6 , 1501)  I TYPE 

FORMAT ( IX, ' NOT  IMPLEMENTED  IN  NORMAL.  I TYPE  =',I5) 

STOP 

ELL 

DX  =  0.  DO 
AN  =  0. DO 

C2  =  FPD(K+8)**2/4.  DO 
DO  1620  1=1,3 
IA  =  I  +  K  +  1 

X( I)  =  XP(I)  -  . 5D0*(FPD(IAj  +  FPDCIA+3)) 

H( I )  =  .  5D0*(FPD(IA+3)  -  FPD(IA)) 

DX  =  DX  +  X( I)*H( I ) 

CONTINUE 
DO  1630  1=1,3 

UN( I )  =  C2*X(I)  -  DX*H(I) 

AN  =  AN  +  UN( I )*UN( I ) 

CONTINUE 

IF( AN  .GE.  l.D-12)  THEN 
AN  =  DSQRT(AN) 

DO  1635  1=1,3 

UN( I )  =  UN( I )  /AN 
CONTINUE 
GO  TO  2000 
ENDIF 

WRITE( 6,1651)  NBO , LSURF , ( XP( I ) , 1=1 , 3 ) 

FORMAT( IX, ' ROUND-OFF  ERROR  IN  NORMAL.  NBO=’ , 15 , ’ LSURF=' , 15 , 
’XP=' ,3E12.  5) 

STOP 

BOX 

L  =  1  +  ( LSUR  -  1 ) / 2 

11  =  IBOXl(L) 

12  =  IBOX2(L) 

DO  1720  1=1,3 

X(I)  =  FPDCK+Il+I) 

H ( I )  =  FPDCK+I2+I) 

CONTINUE 
AN  =  0.  DO 
DO  1730  1=1,3 

J  =  MOD( 1,3)  +  1 

JJ  =  MOD( 1+1,3)  +  1 

UN( I )  =  X( J)*H( JJ)  -  X( JJ)*H( J) 

AN  =  AN  +  UN( I )*UN( I ) 

CONTINUE 

IFCAN  .LT.  1.0D-12)  GO  TO  1650 
AN  =  DSQRT(AN) 

DO  1735  1=1,3 

UN( I )  =  UN( I )/AN 
CONTINUE 
GO  TO  2000 
RPP 

DO  1905  1=1,3 
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UN( I)  =  0.  DO 
1905  CONTINUE 

GO  TO  (1910,1910,1920,1920,1930,1930),  LSUR 
1910  UN( 2)  =  1. DO 
GO  TO  2000 
1920  UN( 1 )  =  1.  DO 
GO  TO  2000 
1930  UN( 3)  =  1. DO 
2000  XMU  =  0.  DO 

DO  2005  1=1,3 

XMU  =  XMU  +  WT(I)*UN(I) 

2005  CONTINUE 

XMU  =  - 1. DO*XMU 
C  AN  =  SIGN( 1.  DO ,  XMU) 

IF(XMU  .LT.  0.  DO)  THEN 
AN  =  -1.  DO 
ELSE 

AN  =  1.  DO 
END  IF 

XMU  =  - 1.  DO*XMU 
DO  2010  1=1,3 

UN ( I )  =  AN*UN ( I ) 

2010  CONTINUE 
RETURN 
END 

FUNCTION  IBFIND(X,MX,XB) 

C  ***  IBFIND  FINDS  CORRECT  TALLY  BINS  *** 

IMPLICIT  REAL* 8  (A-H.O-Z) 

DIMENSION  XB(IOO) 

IBFIND=0 

IF(X  . LT.  XB( 1) )  GO  TO  7777 
DO  1000  1=1, MX 

IF( X  . LE.  XB( 1+1) )  GO  TO  1010 
1000  CONTINUE 

GO  TO  7777 
1010  IBFIND=I 
7777  RETURN- 
END 

SUBROUTINE  CKVCDF( BETA ,NMTRLS , VARIND , WVLNTH , CDF) 

C 

C  THIS  SUBROUTINE  CREATES  A  TABLE  OF  INDEX  OF  REFRACTION  VS 
C  FREQUENCY  VALUES  AND  BUILDS  THE  CDF  FOR  CERENKOV  PRODUCTION 
C  BY  WAVELENGTH. 

C 

IMPLICIT  REAL*8  (A-H.O-Z) 

INTEGER  I , N , MTRL , NMTRLS 

DIMENSION  WVLNTH ( 100) ,CDF(20 , 100) ,VARIND( 100 , 100) , PDF( 100) 

C  ASSIGN  VALUES  TO  CONSTANTS 
N  =  27 

PI  =  DACOS( -1. DO) 

ALPHA  =  1. DO/137. DO 

C  CALCULATE  THE  NUMBER  OF  CERENKOV  PHOTONS  PER  PATH  LENGTH  PRODUCED 
C  FOR  EACH  WAVELENGTH  INTERVAL.  ALSO  CALCULATE  THE  TOTAL  NUMBER 

C  OF  PHOTONS  PRODUCED  PER  PATH  LENGTH 

DO  450  MTRL=1, NMTRLS 
PTGT  =  0.  DO 
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DO  200  1=2, N 

PDF(I-l)  =  2.  DO*PI*ALPHA*DABS( ( 1. D0/WVLNTH( I - 1) 

1  -  1.  DO/WVLNTH(I))* 

2  (1.  DO-1.  DO/BETA/BETA/ VARIND(MTRL, I - 1 ) /VARIND( MTHL ,  I  - 1 )  )  ) 
PTOT  =  PTOT  +  PDF(I-l) 

200  CONTINUE 

C  NORMALIZE  PDF 

DO  300  1=1, N-l 

PDF(I)  =  PDF( I ) /PTOT 
300  CONTINUE 

C  CREATE  CDF  TABLE 

CDF(MTRL, 1)  =  PDF(l) 

DO  350  1=2, N-l 

CDF(MTRL, I )  =  CDFIMTRL, 1-1)  +  PDF(I) 

350  CONTINUE 

450  CONTINUE 
RETURN 
END 

SUBROUTINE  CKVWVL( CDF ,MTRL,NCBIN) 

C 

C  THIS  SUBROUTINE  SAMPLES  THE  WAVELENGTH  OF  THE  EMITTED  CERENKOV 
C  PHOTON  AND  RETURNS  THE  BIN  NUMBER  JN  NCBIN.  IT  REQUIRES  INPUT 
C  OF  THE  ARRAYS  WVLNTH  AND  CDF  WHICH  MUST  BE  OBTAINED 
C  BY  A  PREVIOUS  CALL  TO  CKVCDF. 

C 

IMPLICIT  REAL* 8  (A-H.O-Z) 

INTEGER  I , NCBIN, NBINS,MTRL 
DIMENSION  CDF( 20 , 100) 

DOUBLE  PRECISION  DSEED,  IRAN 
COMMON  /VAXRAN/  IRAN 
C  GET  RANDOM  VARIABLE 

RANF(.  DSEED)  =  GGUBFS(DSEED) 

RND 1  =  RANF( IRAN) 

C  RANDOMLY  SAMPLE  THE  CDF  FOR  CERENKOV  RADIATION  AND  DETERMINE 
C  THE  BIN  NUMBER  OF  THE  SAMPLED  CERENKOV  PHOTON 
NS  INS  =  27-1 
DO  100  1=1 ,NBINS 

IF(RND1  .  LE.  CDF(MTRL,I))  THEN 
NCBIN  =  I 
GO  TO  200 
END  IF 
100  CONTINUE 
200  RETURN 
END 
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APPENDIX  D.  SAMPLE  CROSS  SECTION  INPUT  FILE 

ENERGY  17. 

MATERIAL  AL 
MATERIAL  FE 

MATERIAL  C  . 2729  0  .  7271 
GAS 

DENSITY  . 003986 

MATERIAL  SI  .4675  0  .5325  DENSITY  2.5 

MATERIAL  PB 

TITLE 

. ITS(KEY, ACCEPT)  TEST  CERENKOV  89 . 
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APPENDIX  E.  SAMPLE  CERENKOV  INPUT  FILE  WITH  GEOMETRY 

AND  ZONE  DESCRIPTIONS 


ECHO  1 
TITLE 

TEST  CERENKOV  89 
ELECTRONS 
ENERGY  16.  7 
POSITION  0.  0  0.  0  0.  0 
DIRECTION  0.0  0.0 
HISTORIES  10000 
CUTOFFS  .06819  .06819 
ELECTRON -ESCAPE 
NBINE  10 
NBINT  10 
PHOTON -ESCAPE 
NBINE  10 
NBINT  10 
ELECTRON -FLUX  9  18  NBINE  20 


PHOTON -FLUX  9  18 

NBINE  20 

GEOMETRY 

0.  0 

1.0 

RCC 

0.  0 

0.  0 

-1.  0 

0.  0 

RCC 

6.  35 

0.  0 

0.  0 

0.  0 

0.  0 

0.  0 

0. 0508 

RCC 

3.  81 

0.  0 

0.  0 

0.0508 

0.  0 

0.  0 

0.  635 

RCC 

3.  81 

0.  0 

0.  0 

0. 6858 

0.  0 

0.  0 

0. 00762 

3.  81 

51.  11 

RCC 

0.  0 

0.  0 

0.  0 

0.  0 

0.  0 

6.  35 

51.  11 

RCC 

0.  0 

0.  0 

0.  0 

0.  0 

0.  0 

RCC 

6. 1134 
0.  0 

6.  0 

0.  0 

38. 01842 

-0. 00113 

0.  0 

0.  0 

0. 001130 

3.  0 

RCC 

0.  0 

6.  35 

0.  0 

51.  11 

0.  0 

0.  0 

RCC 

0.  0 

6. 1134 

0.  0 

38.  01842 

36.  744 

0.  0 

0.  0 

0.  0 

RCC 

0.  0 

0.  0 

38. 01842 

36.  744 

6.  35 

13. 785 

RCC 

26. 584 

0.  0 

38. 01842 

0.  0 

0.  0 

6.  1134 

13. 785 

RCC 

26. 584 

0,  0 

38. 01842 

0.  0 

0.  0 

6.  35 

SPH 

28. 716 

0.  0 

38. 01842 

9.  49 

0.  0 

0.  0 

RCC 

18. 504 

0.  0 

38. 01842 

3.  0 

RCC 

4.  95 

26. 584 

0.  0 

38. 01842 

0. 00113 

0.  0 

-0.  00113 

RCC 

6.  0 

26. 584 

0.  0 

51.  80342 

0.  0 

0.  0 

0.  318 
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o-F'j>oo4>orooroof'oo.p'-u>roo^-‘0 


END 


Z1 

+1 

Z2 

+2 

Z3 

+3 

Z4 

+4 

Z5 

+6 

-2 

-3 

-4 

Z6 

+7 

Z7 

+8 

Z8 

+5 

-6 

-9 

Z9 

+9 

+5 

-6 

Z10 

+10 

-5 

-9 

-11 

Zll 

+9 

-5 

-14 

-15 

Z12 

+12 

-10 

-1 

-2 

Z13 

+11 

-10 

Z14 

+13 

+14 

Z15 

+  14 

-13 

Z16 

+10 

+11 

-9 

Z17 

+15 

Z18 

+  16 

Z19 

-1 

-5 

-8 

-10 

END 

*  MAT  ECUT  PTCZ 
0  .  1 


12.  0 
12.  0 
12.  0 


*  The  two  values  immediately  following  the  key  word,  CERENKOV,  allow 

*  the  user  to  set  the  pressure  of  air  or  carbon  dioxide  gas 

*  within  the  geometry,  and  the  temperature  in  kelvin. 

*  (Other  gasses  are  not  included  in  the  current  version. ) 

*  The  pressure  in  units  of  standard  atmospheres. 

*  Once  set,  the  pressure  applies  to  all  of  that  type  of  gas 

*  throughout  the  geometry.  The  temperature  applies  to  all  gasses. 

*  The  code  uses  these  values  to  scale  the  indices  of  refraction 

*  for  air  and  carbon  dioxide. 

*  (The  user  is  still  required  to  generate  the  appropriate 

*  cross  section  tape  for  the  gasses  at  these  new  densities. ) 

*  AIRP  C02P  temp  kaltot 

*  (atm)  (atm)  (kelvin) 

CERENKOV  1.00  2.00  293.0  15 

*  The  material  number  below  is  not  the  same  as  the  material  number 
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*  used  above  in  the  MAT  section.  Instead  it  refers  to 

*  materials  for  which  the  Cerenkov 


*  subroutine  has  index  of  refraction  data.  These  are  currently: 

*  1  Vacuum 

*  2  Air 

*  3  Carbon  Dioxide 

*  4  --  Glass 

*  5  Crystal  Quartz 


*  Note  that  Vacuum  is  used  whenever  we  want  to  set  the  index  of 

*  refraction  to  a  value  of  1. 0. 


MATERIAL 

LAMBDA 

MIRROR 

DIFFUSE 

ACCEPT 

1 

0.  0 

0.  0 

0.  0 

1.  0 

1 

1000.  0 

0.  0 

0.  0 

1.  0 

1 

0.  0 

0.  0 

0.  0 

1.  0 

1 

1000.  0 

0.  0 

0.  0 

1.  0 

3 

0.  0 

0.  0 

0.  0 

1.  0 

4 

0.  05 

0.  9 

0.  0 

0.  0 

1 

0.  0 

0.  0 

0.  0 

1.  0 

1 

1000.  0 

0.  0 

1.  0 

0.  0 

1 

0.  0 

0.  0 

0.  0 

1.  0 

1 

1000. 0 

0.  0 

1.  0 

0.  0 

1 

0.  0 

0.  0 

0.  0 

1.  0 

1 

1000.  0 

0.  0 

1.  0 

0.  0 

1 

0.  0 

0.  0 

0.  0 

1.  0 

4 

0.  0 

0.  0 

0.  0 

1.  0 

1 

0.  0 

0.  0 

0.  0 

1.  0 

1 

0.  0 

0.  0 

0.  0 

1.  0 

4 

0.  05 

0.  9 

0.  0 

0.  0 

1 

1000.  0 

0.  0 

0.  0 

1.  0 

1 

0.  0 

0.  0 

0.  0 

1.  0 

*  AL  20MIL  29. 3P  16.  7MEV  400K  HIST  T11X20  KLIDE  15 
A.  PROBLEM  GEOMETRY 

To  describe  the  geometry  of  an  object  in  ACCEPT  the  user  must  execute  four  dis¬ 
tinct  tasks:  [Ref.  11:  p.  73] 

1.  Define  the  location  and  orientation  of  each  solid  geometrical  body  required  for 
specifying  the  input  zones. 


2.  Specify  the  input  zones  as  combinations  of  these  bodies. 

3.  Specify  the  volumes  of  the  input  zones,  if  necessary. 

4.  Specify  the  material  in  each  input  zone. 


The  use  of  combinatorial  geometry  is  based  upon  a  library  of  geometrical  body  types 
from  which  the  user  chooses  a  specific  combination  to  describe  the  body  he  is  modelling. 
These  figures  are  shown  on  the  following  pages.  [Ref.  11:  pp.  74-82] 
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z 


Figure  £-1:  Rectangular  Paralleiepi ped  (RPP) 

b)  Sphere  (S?K)  —  Specify  the  components  of  the  radius  rector  V  tc  the 
center  cf  the  sphere  and  the  radius  R  of  the  spnere.  ~ 


Figure  E-2:  Sphere  (S?H) 

c)  Right  Circular  Cylinder  (RCC)  —  Specify  the  components  of  a  radius 
rector  7  to  the  center  of  one  base,  the  components  of  a  vector  H 
from  the  center  of  that  base  tc  the  center  of  the  other  base,  and 
the  radius  R  of  the  cylinder. 


Figure  £-3:  Right  Circular  Cylinder  (RCC) 


93 


d)  Right  Elliptical  Cylinder  (REC)  —  Specify  the  components  of  a 
radius  vector  7  to  the  center  of  one  of  the  elliptical  bases,  the 
ccaponenta  of  a  vector  H  frca  the  center  of  that  base  to  the  center 
of  the  other  base,  and  The  components  of  two  vectors  R1  and  R2 

defining  the  major  and  minor  axes,  respectively,  of  the  bases.  This 
body  has  not  yet  been  impleaented. 


Figure  E-4:  Right 'Elliptical  Cylinder  (REC) 

e)  Truncated  Right  Angle  Cone  (TRC)  —  Specify  the  components  of  a 
radius  vector  7  to  the  center  of  one  base,  the  components  of  a 
vector  H  frco  the  center  of  that  base  to  the  center  of  the  other 

base,  and  the  radii  R,  and  R„  of  the  first  and  second  bases, 

•  c 

respectively. 


Figure  E~5 :  Truncated  Rignt  Angular  Cone  (TRC) 

f)  Ellipsoid  (ELL)  —  Specify  the  components  of  the  radius  vectors  71 
and  72  to  the  foci  of  the  prolate  ellipsoid  and  the  length  cf  the” 
major- axis  R. 
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Figure  2-3:  Ellipsoid  (ELL) 

g)  Wedge  (WED)  —  Specify  the  components  of  4  radius  vector  y  to  one  of 
the  corners  and  the  components  of  three  mutually  perpendicular 
vectors  a^ ,  »2  and  a^  starting  at  that  corner  and  defining  the  wedge 

such  that- a^- and  a2""are  the  two  legs  of  the  right  triangle  of  the 

wedge. 


Figure  2-7:  Right  Angle  Wedge  (WED) 

h)  Box  (30X)  —  Specify  the  components  of  a  radius  vector  V  to  one  of 
the  corners  and  the  components  of  three  mutually  perpendicular 
vectors  a1 ,  a2  and  a^  starting  at  that  corner  and  defining  a 

rectangular  parallelepiped  of  arbitrary  orientation. 


V 

Figure  £-3:  Box  (BOX) 


95 


i)  Arbitrary  Polyhedron  (AR3)  —  Specify  the  components  of  k  (k  -  5,  7 
op  8)  radius  vectors,  Y1  through  7^,  to  the  corners  of  an  arbitrary 

nonreer.trant  polyhedrorTof  up  to  Tlx  sides,  and  specify  the  Indices 
of  the  corners  of  each  face  by  means  of  a  series  of  four-digit 
floating  point  numbers  between  **i 230 . "  and  "8765."  (enter  zero  fcr 
the  fourth  index  of  a  three-cornered  face).  These  indices  oust 
appear  In  either  clockwise  or  counterclockwise  order. 


Figure  £-9:  Arbitrary  Polyhedron  (AR3) 

E.1.2  Specification  of  Input  Zones  —  Having  defined  the  necessary 
geometrical  bodies,  the  user  must  then  resolve  the  entire  problem  geometry 
into  input  tones  satisfying  the  following  criteria: 


a)  An  input  tone  say  consist  only  of  either  a  single  homogeneous 
material  or  a  void. 

b)  Every  point  of  the  problem  geometry  oust  lie  within  cr.e  ar.d  only 
one  input  sene. 


c)  The  final  input  tone  must  be  a  non-reentrant  vcid  tone  surrounding 
.  the  rest  of  the  problem  geometry;  any  particle  entering  this  tone  is 
tallied  as  an  escape  particle. 


Input  zenes  are  specified  as  appropriate  combinations  of  the  previously 
defined  bodies.  Such  combinations  may  be  as  simple  as  Just  a  single  body, 
or  they  may  consist  of  complex  intersections,  unions  and  differences  of 
various  bodies.  We  illustrate  the  principles  of  input  tone  specification 
with  the  following  examples  where,  for  simplicity,  we  omit  the  escape  tone. 
Each  example  involves  only  two  zones,  A  and  3,  defined  by  the  crcss-r.atchir.g 
in  Fig.  E-1G. 
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INPUT  ZONE  A 


INPUT  ZONE  B 


BODY  #3 


In  Fig.  E-lOa,  zone  A  consists  of  a  sphere,  body  #1  ,  that  is  tangent  to 
zone  3  which  consists  of  a  right  circular  cylinder,  body  92.  Input  zone 
specification  is  simply 


A  -  «-1 

and  3  *  *2  . 

That  is,  incut  zone  A  consists  of  all  soatlal  points  that  lie  within  body 
#1 ,  and  similarly  for  zone  3. 

In  Fig.  E-lOb,  the  sphere  is  inserted  into  a  hole  that  has  been  cut  in 
the  cylinder  so  that 


A  -  *1 

and  3  -  *2  -1  . 

Thus,  incut  zone  3  consists  of  all  soatial  points  that  lie  within  body  92 
AMD  not  within  body  sH  .  Input  zone  3  is  specified  as  the  clfferer.ce  between 
two  do dies. 

In  Fig.  E-TOc,  bodies  #1  and  92  consist  of  the  same  homogeneous 
material  (or  void) ,  but  they  are  imbedded  within  a  second  right  circular 
cylinder,  body  91,  of  another  material.  The  specification  is 

A  -  *1  CR  *2 
and  3  *  *3  “1  -2  . 

Thus,  input  zone  A  consists  of  all  soatlal  points  that  lie  within  EITHER 
body  OR  aody  42.  This  is  an  example  of  input  zone  specification  as  a 
union  of  bodies. 

In  Fig.  E-IOd,  the  intersection  of  body  #1  and  tody  92  consists  of  a 
single  homogeneous  material;  the  rest  of  the  space  within  body  91  is  filled 
with  another  material.  The  specification  is 


A  -  *1  *2 

and  3  -  *3  -1  OR  *3  "2. 

Thus,  ir.cut  zone  A  consists  of  all  soatial  points  that  lie  within  body  9' 
AMD  within  body  92. 


Mote  that; 

a)  The  OR  operator  refers  to  all  following  body  numbers  until  the  next 
OR  operator  is  reached  or  a  new  input  zone  is  initiated. 

b)  The  AMD  operator  is  implied  before  every  body  number  that  is  r.ct 
preceded  by  an  explicit  OR  operator,  except  that  the  firs:  OR 
operator  of  a  union  is  an  implied  EITHER. 

E.1.3  Volume  Specification  —  A  general  scheme  for  the  precise 
calculation  of  the  volumes  of  input  zones  defined  by  the  combinatorial 
method  is  not  possible.  The  user  may  select  one  of  three  ootiens  via 
Cparameter (1 ) ]  associated  with  the  GEOMETRY  keyword.  The  default  value  of  I 
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will  cause  the  code  to  set  all  volumes  to  1.0  ea.^  A  value  of  1  allows  the 
user  to  read  In  the  volumes  as  described  below.  Finally,  If  the  geometry  is 
such  chat  a  satisfactory  method  exists  for  calculating  the  volumes 
Internally,  the  user  may  set  Cparameter( 1 ) ]  equal  to  2  and  use  a  correction 
run  to  insert  the  necessary  logic  at  the  proper  place  near  the  end  of 
Subroutine  JOCZN.  The  proper  location  is  commented  in  the  coding. 

E.1.4  Material  Specification  —  A  material  index  1s  assigned  to  each 
input  zone.  A  zero  index  defines  a  void  zone;  otherwise,  the  material 
indices  are  defined  by  the  order  in  which  the  materials  are  specified  in 
executing  the  cross  section  generating  code  as  described  In  Sec.  2.2.2.  The 
method  of  inputting  the  material  indices  is  described  under  the  keyword 
GtCMETSr. 

E.2  Input  Data 

The  geometry  input  for  the  ACCZPT  eodes  Is  inserted  according  tc  the 
following  sequence.  The  data  is  Inserted  in  free  format  form  with  spaces  or 
commas  as  delimiters.  Simple  examples  can  be  found  In  Secs.  3.3,  3.6  and 
S.8  of  Appendix  3. 

E.2.1  Body  Data  —  The  body  data  begin  immediately  after  the  line 
containing  the  GEOMETRY  keyword.  The  method  of  describing  each  of  the  fcodv 
types  is  discussed  in  Sec.  E.1.1  and  illustrated  in  Table  Z-1 .  The 
description  of  each  new  body  must  begin  a  new  line  of  input,  and  the  first 
parameter  on  that  line  luosc  be  the  appropriate  three  character  code  for  the 
body  type.  Table  E-1  lists  the  additional  input  parameters  required  (no 
defaults)  for  each  body  type  in  their  proper  sequence.  The  user  is  free  to 
distribute  these  parameters  over  am  many  lines  as  he  pleases.  A  line 
beginning  with  the  parameter  END  signals  that  the  description  of  ail  of  the 
problem  bodies  is  complete. 

E«2.2  Input  Zone  Data  —  Geometrical  specification  of  the  input  zones 
begins  immediately  after  the  line  containing  the  E;,D  parameter  for  the  bccy 
data.  The  method  of  describing  the  input  zones  in  terms  of  the  input  bodies 
is  discussed  in  Sec.  E. 1.2.  Body  numbers  are  determined  by  the  order  in 
which  the  bodies  are  read  in.  The  description  of  each  new  input  zone  must 
begin  on  a  new  line  of  input,  and  the  first  parameter  on  that  line  must  be  a 
character  string  beginning  with  the  letter  Z.  This  parameter  is  followed  by 
*  string  of  parameters  that  specifies  the  input  zone  following  the  form  of 
the  right  hand  sides  of  the  equations  of  Sec.  E.1.2.  For  example,  input 
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Table  2-1  Data 

Required 

to  Oescri 

bed  Eacn  Body  Type 

Body  Type 

Real 

Data  Def 

L.ning  Particular  Body 

BOX 

Vx 

Vy 

Vz 

Alx 

Aiy 

A1Z 

A 2x 

A 2-/ 

A2z 

A3x 

A3y 

A3  = 

RPP 

Xmin 

Xmax 

Tain 

Tmax 

Zmin 

2  max 

SPH 

Vx 

Vy 

V: 

?. 

RCC 

Vx 

Vy 

Vz 

Hx 

Hy 

Hz 

R 

re: 

Vx 

Vy 

vz 

Hx 

Hy 

Kz 

R1  x 

Ri  y 

Riz 

R2x 

R2y 

R2z 

r* r 

71  x 

9 

VI  y 

viz 

V2x 

72  y 

72  z 

TRC 

Vx 

Vy 

vz 

Hx 

Hy 

i;  , 

R1 

R2 

WE  3 

Vx 

Vy 

vz 

A1  x 

Aiy 

Alt 

kZx 

A2y 

A2s 

A-x 

A3y 

■A3b 

ARB 

VI X 

VI  y 

Viz 

V2x 

V2y 

V2z 

V3x 

v  3y 

V2z 

VUX 

VUy 

VUz 

V5x 

V5y 

V5  z 

V6x 

V6y 

76  z 

V7x 

V7y 

V7z 

V8x 

VSy 

VS: 

Face 

Oescri pti 

ens  (se 

e  note 

below) 

~10  .Vo  Oata 


Note: 


The  final  Line  of  the  arbitrary  polyhedron  input  contains  a  four- 
Oigit  real  number  for  each  of  the  races.  Thirty  data  values  are 
required  f or  this  body  type;  if  there  are  fewer  than  eight  corners 
and  six  faces,  zero  values  .oust  be  entered.  See  Section  E.’.t  for 
mere  details. 
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lines  describing  the  two  input  zones  in  rig.  E-10d  are 
2001  *1  *2 

2002  *3  “1  OR  *3  '2  . 

The  user  is  free  to  distribute  the  parameters  necessary  for  describing  an 
input  zone  over  as  many  lines  as  he  pleases.  A  line  beginning  with  the 
parameter  END  signals  that  the  description  of  all  of  the  problem  input  zones 
is  complete. 

E.2.3  Volume  Data  —  If  [parameter ( 1 ) ]  associated  with  the  GEOMETRY 
keyword  is  equal  to  1 ,  the  array  containing  the  volume  data  is  inserted 
immediately  after  the  line  containing  the  END  parameter  for  the  data 
specifying  the  geometry  of  the  input  zones.  The  input  zones  are  numbered 
according  to  the  order  in  which  they  are  read  in.  The  array  must  contain  an 
entry  for  each  input  zone  (no  defaults).  If  [parameter ( 1 ) ]  is  not  equal  to 
1  ,  these  input  data  are  omitted. 

E.2.9  Material  Data  —  For  input  of  material  data,  return  to  the 
discussion  under  the  GEOMETRY  keyword. 
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APPENDIX  F.  PARTICLE  WEIGHT 

ITS  uses  a  Monte  Carlo  approach  that  often  does  not  exactly  simulate  physical 
transport.  For  instance,  each  ITS  particle  might  represent  a  number  W  of  particles 
emitted  from  a  source.  Here,  W  is  the  initial  weight  of  the  ITS  particle.  In  an  exact 
simulation,  the  W  physical  particles  would  each  have  a  different  random  walk,  but  the 
ITS  particle  representing  these  W  physical  particles  will  only  have  one  random  walk. 
Many  particles  of  weight  W  are  initiated  and  followed  producing  statistically  identical 
results.  Although  this  is  not  an  exact  simulation,  weight  is  preserved  in  ITS  in  the  sense 
of  statistical  averages  and  therefore  in  the  limit  of  large  particle  numbers  (including 
particle  production  or  losses  that  may  occur).  Each  ITS  particle  weight  is  tallied  so  that 
the  final  results  (tallies)  reflect  the  results  expected  from  W  physical  particles.  (Weights 
do  not  have  to  be  integers!)  This  allows  users  to  normalize  the  output  to  the  desired 
source  strength.  ITS  results  are  normalized  to  the  number  of  incident  source  particles. 
[Ref.  S:  pP.  53-54] 

As  an  example,  the  Cerenkov  routines  in  this  work  calculate  the  total  number  of 
Cerenkov  photons  produced  for  each  electron  substep  (small  pathlength).  Rather  than 
tracking  each  of  these  photons  individually,  a  random  position  along  the  substep  is  se¬ 
lected  as  a  representative  production  site  and  a  Cerenkov  photon  of  appropriate  weight 
is  generated.  This  photon  weight  represents  all  Cerenkov  photons  generated  over  this 
pathlength.  This  weight  is  then  multiplied  by  the  parent  electron  weight  (corresponding 
to  the  number  of  electrons  represented  by  the  ITS  electron  taking  the  substep),  and  the 
final  weight  is  assigned  to  the  Cerenkov  photon.  The  photon  is  then  ray  traced  through 
the  problem  geometry.  Then  the  next  electron  substep  is  taken  and  the  process  is  re¬ 
peated. 

The  utility  of  particle  weight  does  not  rest  solely  upon  normalizing  the  source. 
Weight  is  a  critical  tooi  used  in  biasing  problems  to  improve  computational  efficiency. 
For  example,  a  user  may  be  interested  in  studying  the  results  of  an  interaction  that  has 
a  small  probability  to  occur.  Normally,  this  interaction  is  seldom  sampled  and  necessi¬ 
tates  large  numbers  of  histories  to  get  a  statistically  significant  result.  A  better  approach 
would  be  to  increase  the  sampling  frequency  for  longer  wavelengths  with  a  correspond¬ 
ing  weight  adjustment  in  the  interaction.  This  can  be  done  efficiently,  providing  a  higher 
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figure  of  merit,  if  care  is  taken  to  alter  the  results  so  that  the  biasing  is  removed  before 
scoring. 

Summarizing,  weight  represents  the  ITS  particle's  relative  contribution  to  the  final 
tallies.  Its  magnitude  ensures  that  the  expected  physical  results  are  preserved  in  the 
sense  of  statistical  averages  and  therefore  in  the  limit  of  large  numbers  of  ITS  particles. 
Weight  enables  the  user  to  control  the  number  of  particles,  sampling  the  critical  part  of 
a  problem  to  improve  the  accuracy  of  selected  results.  Biasing  used  with  weight  control 
can  often  yield  significantly  improved  figures  of  merit.  [Ref.  8:  p.  54] 
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